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ABSTRACT 


This  report  contains  a  description  of  the  models  to  be  used  in 
analyzing  the  capabilities  of  ground-based  sensors  in  determining  the 
mass  of  orbiting  bodies,  a  description  of  the  model  coefficients,  and  the 
justification  for  their  selection.  Relations  are  derived  for  computing 
sensitivity  coefficients  and  their  coupling  to  mass  variance. 


SECTION  I 


INTRODUCTION 


This  document  is  Technical  Report  Number  1,  the  first  in  a  series  of 
three  theoretical  reports  prepared  under  Air  Force  Contract  F196 28-67-00041. 
It  is  concerned  with  the  mathematical  models  and  relationships  necessary 
to  perform  a  detailed  maximum-likelihood/minimum-variance  error  analysis 
of  the  capability  of  ground-based  sensors  in  determining  the  mass  of  a 
satellite.  It  is  specialized  to  account  for  the  following  restrictions: 

1)  The  sensors  observe  the  satellite  without  error. 

2)  The  satellite  is  a  sphere  of  5-meter  diameter. 

3)  All  the  physical  characteristics  of  the  satellite  except  mass 
are  known  without  error. 

4)  All  the  error  in  the  computer  mass  results  from  errors  and 
uncertainties  in  the  knowledge  of  the  orbit-perturbing  forces. 

Technical  Report  Number  3  will  be  a  companion  document,  extending  the 

theoretical  development  to  remove  the  restrictions  of  perfect  sensor 
observations  and  perfect  knowledge  of  the  non-mass  body  characteristics. 

In  addition,  the  body  shapes  will  be  generalized  from  the  sphere  of  this 
report  to  include  also  one  of  a  pair  of  tumbling  cylindrical  objects  of 
length  10  meters  and  diameters  2  and  5  meters,  respectively. 

The  theoretical  basis  for  mass  determination  rests  in  the  fact  that  the 
mass  of  a  man-made  satellite  drops  out  of  the  gravitationally  induced 
motion  of  the  satellite  and  appears  only  in  the  effects  produced  by  the 
non-conservative  force  fields.  For  the  sort  of  satellites  specified  in 
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the  contract ,  solar  radiation  pressure  and  atmospheric  drag  are  the  only 
two  non-conservative  phenomena  whose  effects  are  reasonably  observable. 
Hence  the  problem  is  to  separate  the  gravitationally  induced  motion  from 
the  net  motion  and  then  use  the  proper  solar  pressure  and  atmospheric 
drag  models  to  extract  mass-parameter  values  for  the  satellites  tracked. 

This  contract  is  not  concerned  with  the  processing  of  real  tracking 
data.  Rather  it  is  a  study  to  ascertain  how  accurately  mass  parameters 
can  be  typically  determined  in  the  way  just  described  and  what  are  the 
critical  error  sources  in  such  a  determination. 

In  choosing  the  mathematical  models  upon  which  to  base  the  necessary 
orbital  calculations,  one  finds  oneself  deeply  involved  with  questions  of 
practical  computation.  A  central-force  gravity  law,  for  example,  leads 
to  expressions  for  the  state  of  the  orbit  which  are  closed-form  functions 
of  eccentric  anomaly  and  hence  are  easily  computable.  Unfortunately, 
this  solution  is  far  too  inaccurate  a  representation  of  the  real  world 
to  be  used  as  is  in  orbit  or  vehicle-parameter  determination. 

More  exact  models  of  satellite  dynamics  do  not  lead  to  closed- form 
solutions.  Not  only  must  the  trajectories  that  evolve  from  the  direct 
use  of  these  more  accurate  models  be  determined  by  numerical  integration, 
but  also  the  state  and  mass-parameter  estimates  and  the  associated  sensi¬ 
tivity  matrices  must  be  determined  by  numerical  integration. 

To  handle  this  problem,  the  use  and  extensions  of  techniques  contained 
in  NASA's  MINIVAR  family  of  orbit-determination  computer  programs  will  be 
made  in  this  contract. 
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In  what  follows  in  this  first  Technical  Report,  the  basic  theory  of 
two-body  mechanics  will  be  presented,  the  means  for  correcting  this  model 


for  various 
the  effects 
as  to  avoid 
selves  will 


pertubation  effects  will  be  discussed,  and  methods  for  computing 
of  variations  in  the  model  coefficients  will  be  developed  so 
the  need  for  numerous  integrations.  Finally,  the  models  them- 
be  discussed. 
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SECTION  II 


TWO-BODY  MECHANICS 


A  two-body  or  Keplerian  earth  orbit  leads  to  an  essentially  closed- 
form  expression  for  the  state  of  the  orbit  as  a  function  of  the  state 
at  any  prior  time . 

In  particular,  if 


r 


x 


y 


_z . 

is  the  column  vector  of  satellite  position  in  the  geocentric  inertial 
coordinate  system  (X,  Y,  Z)  depicted  in  Figure  1,  then  the  central- 
force-law  acceleration  is 


d^r/dt2  =  -  /«r/r\  (2.1) 

where  ju  is  the  gravitational  constant  of  the  earth.  Given 
initial  conditions  at  some  time  tQ  for  position  and  velocity,  r^  and 


North 

Pole 
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mean  vernal  equinoxy 
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X 
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Equatorial 
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0“  January  1 


Figure  1.  Geocentric  Inertial  Coordinates 
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fo  =  (dr/dt)Q,  respectively,  one  can  write  the  position  and  velocity 
at  any  later  time  t  as  ^  3>  4* 


where 


and 


r  =  rQ  f(4E)  +  fQ  g(4E), 

dr/dt  =  rQ  fj.(A  E)  +  tq  g^AE),  (2.2) 

f(AE)  =1  -  (a/rQ)(l  -  cos  A  E), 

g(AE)  =  vQ{a.//<  sin  AE  +  (ado//<)(l  -  cos  A  E), 

ft(AE)  =  -  (a/r)(l/ro)(/</a)^1//2^  sin  A  E, 
gt(A  E)  =  1  -  (a/r)(l  -  cos  A  E), 

a  =  ^rQ/(2//  -rQr  ),  the  semi -major  axis 

d  =  r  *r  , 
o  -o  ~o’ 

AE  =  E  -  Eq,  the  change  in  the  eccentric  anomaly  from  tine  t  , 
r/a  =  (1  -  cos  A  E)  +  (rQ/a)  cos  AE  +  dQ(/<  sin  AE. 


Kepler's  equation  provides  the  means  for  finding  A  E,  given  EQ,the  time  step 
at,  and  the  eccentricity  e: 

n  A  t  =  AE-e  (sin  E  -  sin  Eq), 
or,  after  some  elementary  substitutions, 

nAt  =  AE  +  (do/x/^Ta)(l  -  cos  A  E) 

+  (rQ/a  -  1)  sin  AE,  (2.4) 

where 

At-  =  t  -  t 

o 

n  =  the  mean  motion,  ( ju  . 

^Superscript  numerals  denote  entries  in  the  References  Section  of  this  report. 
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The  solving  of  (2.4)  for  AE  completes  the  two-body  solution. 


For  notational  simplicity,  the  equations  (2.2)  can  be  expressed 
in  state-variable  form: 

X  =  x(xq,  A  E) ,  (2.5) 

where 


"1 


x  •  I  • 

-  I  • 


;  xx  =  x,  x2  =  y,  x3  -  z,  x4  -  x,  x5  =  y,  x6  -  z, 


Lx6 


X  = 

-0 


Xol] 


X-/ 

L  06J 


;  x^-j  -  x  ,  etc, 
9  Ol  o> 


Hence  a  small  change  &  x  in  the  state  caused  by  a  variational  change 
Ax^  in  the  initial  conditions  is,  to  first-order  accuracy, 

Ax  =  $(t,  to)  Axo,  (2.6) 

where 

|(t,  to)=  [dx./dxo.  _  . 

The  36  elements  of  the  transition  matrix  l(t,  t  )  are  derived  in  Appendix 
IV.  They  are  evaluated  on  the  nominal  two-body  orbit  calculated  from  (2.2). 


-  6  - 


SECTION  III 


PERTURBATIONS  TO  TWO-BODY  MOTION 


Although  the  means  has  just  been  developed  for  accounting  for 
variations  in  initial  conditions,  two-body  assumptions  are  still  inadequate. 
The  closed-form  solutions  presented  above  will,  however,  be  useful  in 
the  more  sophisticated  computations  which  can  be  made. 

A .  GENERAL  EQUATIONS  OF  MOTION 
In  general, 


dx/dt  =  f(x,  t)  +  F(x,  £,  u,  t), 
djf/dt  =  G(x,£  ,  u,  t), 


(3-D 


(3.2) 


where  x  is  as  defined  for  equation  (2.5),  5  Is  a  six-dimensional  state 
vector  describing  rigid-body  rotation,  having  the  Euler  angles  of  the 
satellite  as  its  first  three  elements  and  their  time  derivatives  as  its 
remaining  three,  and  u  is  a  vector  containing  all  the  model  parameters 
about  which  there  exists  modeled  uncertainty. 

The  vector  f(x,  t)  describes  the  two-body  accelerations,  while  the 
non-Keplerian  perturbations  enter  with  F(x,£ }  u,  t),  which  contains  also 
all  the  dynamic  biases  and  uncertainties: 


(3-3) 


The  rotational  states  £  arise  for  aspherical  satellites  to  describe  the 
cross  sectional  area  presented  to  the  drag  media  and  to  account  for  the 
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dynamic  coupling  between  rotational  and  translational  energy. 


B.  EQUATIONS  FOR  SPHERICAL  SATELLITES 

For  the  purposes  of  this  Task  I  effort,  only  spherical  satellites 
are  to  be  considered.  To  a  high  degree  of  accuracy,  then, ^  may  be 
dropped  from  (3»l)  and  ( 3 -3 ) ?  and  equation  (3-2)  can  be  discarded. 

Only  the  following  need  therefore  be  considered: 

dx/dt  =  f(x,  t)  +  F(x,  u,  t),  (3*4) 

The  non-Keplerian  perturbations  contained  in  F(x,  u,  t)  will  consist  of 
l)  an  atmospheric  drag  acceleration 

”4rag  "fV1)  +  u2^  (3-5) 


where 


u,.  =  (A/m)C^,  the  ballistic  coefficient,  constant  but  unknown, 
u^(t)  =  a  stationary  random  error  in  the  ballistic  coefficient  with 

2  —  It I 

autocorrelation  function  <j,  e  1  /tai  ,  due  to 

drag  * 

atmospheric  density  uncertainties, 


g.  (x,  t) 


3J 


=  -  | +  uex2)2  +  (x5  -  +  xe 


X. -Ho  x„ 

4  e  2 
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where,  in  turn, 


^  =  mean  atmospheric  density  at  x  and  t, 

a)  =  rotation  rate  of  the  earth: 
e  1 

2)  a  solar  pressure  acceleration  (see  Appendix  I) 


^solar  =  **3 


(3.6) 


where,  for  the  spherical  satellites  in  this  Task  I  effort, 

u^  =  (l+Ak^/9)  ( ^non/cHA/m) ,  a  solar  "ballistic"  coefficient,  constant 
but  unknown, 


and 


Y(x,  t) 


(Res/res)2(Plir+P2i0s) 


> 


where 


I 

nom 

R 

es 

r 

es 


c 


A 


A 

1 

SS 


diffuse  reflectivity  of  satellite, 

solar  irradiance  at  nominal  earth-sun  distance  R  , 

es , 

nominal  earth-sun  distance, 
actual  earth-sun  distance, 

speed  of  light 

unit  vector  from  center  of  earth  to  center  of  satellite, 
unit  vector  from  center  of  sun  to  satellite, 
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f  Xq  (a.  +  a_J  cos  a 


i 


\q  a 


0 


lot!  ^  tt/2  -  B. 

1  '  A 

rr/2  -  ^  |  a  |  ^  tt/2  + 

tt/  2  +  ^  |  a  |  <  TT 


r  (l  ■+  4qaJ 

(1  ,  xq%s) 

0 


1  a  j  £  tt/2  -  B. 


tt/2  -  B^  ^ 


tt/2  4  B 


tt/2  +  B,  \  a|  £ 


TT 


A 

-(i 


cos  a 

sin  a  =  (l-cos^  a  ) 


i  ).  0  £  a  <u  n. 

55  (1/2) 


q  =  earth  albedo , 


R  - 


a-,  = 


ss 

a 


ratio  of  earth  radius  to  distance  between  satellite  and 
earth  center , 

-1, 

COS  A  , 

-  (.0417  +  .5431?0/3 

[.0444  -  3 .17(7'  -.77)3  +  .0045  (X- .77)  sin  [  14  •  3(^-.77)tt])/3 
a2  [  14  s  -  sest}r  -  e’t  (2  +  sy)]  /2 
[  ass+al/2  0  +  1  -  s(l  +  sy)d  ]}  cos  a 


\2  l  ( lA-sin  a)3  +  (\-sin  a)'  J 


(1-X2)3/2 


sm  a 


( l+X2-  27.  sin  a)3/^ 


T  =  -4  4-  9. 34 

y  =  (a-Tr/2)/B? 

s  4’1  ’  yi°’ 

i  i  5  y-o, 

d  =  3-7  +  59(4-. 77)"; 
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3)  an  acceleration  due  to  solar  and  lunar  attraction 


“  A 

i 

se 

A  “I 

i 

83 

+/^hioon 

"a 

i 

me 

A 

^ms 

I*  05 

-sun,  moon  ' sun 

r2 

r2 

R2 

L  em 

"  R2 
ms 

L  es 

ss  - 

(3.7) 


where 


/^sun’/^aoon 


i  ,  $ 
se7  me 


R  ,  R 
es7  em 


A  /\ 

i  ,  i 
ss 7  ms 


R  ,  R 
ss7  ms 


=  gravitational  constants  for  sun  and  moon, 

=  sun-to-earth  and  moon-to-earth  unit  vectors, 

=  sun-to-earth  and  moon-to-earth  distances, 

=  sun-to-satellite  and  moon-to-satellite  unit  vectors, 
=  sun-to-satellite  and  moon-to-satellite  distances; 


k)  and,  finally,  a  geopotential  acceleration  due  to  the  oblateness  (more 
generally,  the  asphericity)  of  the  earth  (see  Appendix  III): 


N 


n 


-(£)  ‘  T  V  Jn  “n”  +  S"“8rad  VnB)> 

n=2  m=0 


(3.8) 


where  T  =  the  rotation  transformation  that  takes  the  earth- fixed 

geocentric  coordinate  system  (X^,  Y^,  Z'*'),  defined  as  the 
right-hand  system  with  X^  at  Greenwich  and  Z^  being  the 
north-directed  polar  axis  at  time  t,  into  the  inertial 
geocentric  coordinate  system  (X,  Y,  Z)  defined  in  Figure 
1, 
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C  ,  S  =  the  geopotential  coefficients  from  (n,  m)  equal  to 
nm  ’  nm 

(2,  2)  upto  (N,  N)  whose  published  values  have  constant 
but  unknown  biases  on  them, 


grad  = 


column  (— y, 


dx 


n 


(/•Vr)  (R/r)n  (sin  3) 


\ 

cos  m X 

4  r 

sin  mX 

V  ) 


where 


M  = 


r  = 


R  = 

F^Csin  3 )  = 
n 

3,  X  = 


earth’s  gravitational  constant , 

the  geocentric  distance  to  the  satellite  at  time  t, 
earth's  mean  equatorial  radius, 
associated  Legendre  polynomials , 

satellite  latitude  and  longitude,  respectively,  at  time  t. 


Hence 

(3.9) 

C .  ENCKE  INTEGRATION 

Given  the  initial  condition  x(t  )  and  the  actual  values  of  the 

—  o 

parameters  that  appear  in  (3*9),  equation  (3-4)  can  be  integrated  numerically 
to  yield  the  satellite  trajectory  x(t).  Since  in  the  class  of  orbits 


F(x,  u,  t)  = 


0 

+ 

0 

0 

tlrag 

TI 

-solar 

p 

-sun,  moon 

+ 


0 

F" 

—0 
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germane  to  this  contract  the  earth's  central-force  field  is  the  strongly 
dominant  effect,  an  Encke  integration  of  (3.4)  will  generally  be  the 
most  accurate  approach  for  a  given  integration  step  size  and  arithmetic 
precision.’*' 

The  Encke  method  that  will  be  used  computes  small  perturbations 
about  the  two-body  orbit  defined  by  equations  (2.2)  and  (2.5).  Call  the 
two  body  orbit  x°(t),  or,  more  exactly, 

x°(t)  =  x°[t,  x(tQ)]  . 

It  satisfies  the  equation 

dx°/dt  =  f(x°,  t),  (3-10) 

x(tQ)  =  given. 

If  we  define 

Ax(t)  =  Ax  [t,  x(tQ)]  =  x(t)  -  x°  [t,  x( tQ )] 

and  subtract  (3.10)  from  (3*4),  we  find 

d(Ax)/dt  =  f(x,  t)  -  f(x°,  t)  +  F(x,  u,  t), 

Ax(tQ)  =  0.  (3.11) 

The  integration  of  this  equation  by  an  appropriate  starting  technique,  such 
as  the  Range-Kutta-Gill,  and  an  appropriate  long-term  technique,  such  as 
an  Adams  interpolation  method,  completes  the  Encke  orbit  computation. 
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Reference  4  presents  the  Runge-Kutta-Gill  and  Adams  methods  that  will 
be  used,  except  that  it  associates  the  Adams  method  with  a  Cowell  integrator. 
The  CPCEI  Detail  Specifications,  Part  II,  for  the  present  contract  will 
describe  the  same  Adams  method  in  the  context  used  here. 

D.  RECTIFICATION 

If  we  continue  assuming  that  x(tQ)  anc*  the  actual  parameter  values 
are  known  exactly,  care  must  nontheless  be  taken  that  the  magnitude  of 
Ax(t)  does  not  exceed  certain  assignable  bounds,  or  the  Encke  integrator 
will  lose  its  accuracy.  If  it  does,  say  at  time  t  ,  then  the  integration 
must  be  stopped  and  the  generated  value 

xUr)  =  X°  [v  *(tQ)]  +  A*  [tr,  x(to)] 

used  as  a  new  initial  condition.  A  new  two-body  trajectory  is  generated 
from  (2.2)  and  (2-5), 

x°(t)  =  x°  [t,  x(tr)]  > 

a  new  Ax(t)  defined, 

Ax(t)  =  AX  [t,  x(tr)]  =  x(t)  -  x°  j_t,  x(tr)]  , 

and  the  process  continued.  The  technique  is  known  as  rectification,  and 

proceeds  in  one  dimension  as  shown  in  Figure  2,  below.  The  times 

t  ,  t  .....  are  the  rectification  times, 
r  *  r  9  y 
1  2 

The  generation  of  an  orbit  under  the  assumption  of  perfect  knowledge 

O  1 

is  carried  out  in  the  Reference  Mode  of  MINIVAIv  ’  in  the  manner  discussed 
in  this  section. 
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Figure  2.  Actual  Orbit  as  It  Would  Be  Represented  in  an  Encke  Integration 


SECTION  IV 


ESTIMATION  AND  PREDICTION 


The  assumptions  of  the  previous  section,  i.e.,  that  the  initial 
conditions  and  the  model  parameters  are  completely  known,  led  to  an 
Encke  integration  which,  in  reality,  only  God  can  perform.  Rather,  a 
ground-based  observer  has  incomplete  knowledge  of  the  initial  conditions 
and  model  parameters  and  often  can  see  only  some  nonlinear  combination  of 
some  of  the  satellite  states.  Even  that  which  can  be  seen  is  masked  in 
uncertainty  because  of  fluctuating  measurement  noise  and  fixed,  but 


measurement  biases. 


A.  FR0BLEI1  FORMULATION 


Consider  the  equations 


d(Ax)/dt  =  f  (x,  t)  -  f(x°,  t)  +  F(x,  u,  t), 


(4.1) 


Az  =  h(x,2|  ,  b  ,  t)  -  h  (x°,  0,  0,  t), 


(4.2) 


where  Ax(t)  is  the  Encke  variation  from  the  two-body  orbit  x°(t)  discussed 
in  Section  III^  and  the  vector 


z  =  h(x>  b  ,  t) 

defines  the  observations  .z(t)  on  the  orbit,  where 


(4.3) 


23  (t)  =  zero-mean  Gaussian  white  noise, 


b  =  constant  but  unknown  biases , 

A^(t)  =  z  -  z°,  the  difference  between  the  actual  noisy  observation 


and  the  noise-free  observation  of  the  two-body  motion. 
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For  the  special  case  of  perfect  sensors  to  be  considered  in  this  Task 
I  report,  z(t)  is  always  h(x,  0,  0,  t),  which  we  will  shorten  to 

z  =  h(x) .  (4.3  ) 

In  addition  to  the  physically  meaningful  equations  (4.1)  and  (4*2), 
additional  equations  can  be  written  which  describe  how  the  unknown 
parameters  and  biases  within  F(x,  u,  t)  change  with  time.  Three  of  these 
will  be  modeled  and  estimated  from  received  sensor  data:  the  three  "ballistic" 
coefficients  u^,  vl^,  and  u^  in  expressions  (3*5)  and  (3*6),  which  have  the 
dynamics 

du^/dt  =  -(l/rd)  \  +  w^t), 
du^dt  =  0 , 

dv^/dt  =  0,  (4.4) 

where  w^(t)  is  a  zero-mean,  Gaussian  white-noise  process  with  power 

O 

^^Td^drag  Per-1^t-d°uble-bandwidth  in  rad. /sec. j  i.e.,  it  has 
covariance 

Cov  [wx(t),  wl(t’)]  =  ( 2/ 'fc’d)  <T drag  <^ ( t-t ' ) ,  (4.5) 

where  </(t)  is  the  Dirac  delta  function >  and  where  in  general,  for  random 
vectors  a(t),  b(t)^ 

Cov  [a(t),  b(tT|  =  Efa(t)  b(t)]  -  E[a(t)[ 
where  E(  )  denotes  mathematical  expectation  and  (  )  denotes  matrix 

transpose. 


--IT-- 


Additional  uncertain  parameters  are  contained  in  the  oblateness 
contribution  F  to  F(x,  u,  t).  These  are  the  geopotential  coefficients 


J  =  column  (C20,  C^Q, . . . 


which  are  constant  but  have  unknown  biases  on  their  published  values.  These 
errors  will  be  modeled,  but  sensor  data  will  not  be  used  to  correct  the 
published  values . 

For  simplicity  of  notation,  the  vector  v^  will  be  used  to  denote 
those  biases  and  parameters  which  will  be  actively  estimated,  and  v  will 
be  used  to  account  for  those  that  will  be  modeled,  but  not  actively 
estimated.  If  y  is  defined  as  the  vector  which  contains  all  variables  that 
are  to  be  actively  estimated,  then  it  is  the  9-vector 


Lx 


Z  = 


*1 


(4.6) 


where 


Now,  for  any  random  vectors  a(t)  and  b(t)  we  will  define 

§(n/k)  =  a(t  /t,  ),  the  minimum  -  variance  estimate  of  the  vector 
n  k 

a  at  time  t  =  t  based  upon  data  upto  and  including  time 


n 


t  =  V 


Pab(n/k)  =  CovQi(n)  ~'a(n/k),  b(n)  -  'b(n/k)|  ,  the  covariance  matrix 
on  the  estimation  errors  between  a  and  b, 
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and  the  short- form  notation. 


Cov 


£  a(t),  a(t)] 


In  this  notation,  the  problem  of  minimum  -  variance  estimation  on  the 
vector  y  is  simply  the  problem  of  choosing  a  filter  for  the  data  sequence 
^Az(n)}  such  that 

L  a (n/n)  =  trace  i AP  (n/n)A*\  =  min.,  (4.7) 

t  yy  j 


where* 


<T  t  (n/n) 

•'l 


a. 

1 


“  Cov  Jy^(n)  -  y^(n/njj 


=  a  pre-as signed  weight  given  to  an  error  in  the 
estimate  of  the  ith  variable  y.  , 


“l 


°9 


trace  (  )  =  the  sum  of  the  diagonal  elements  of  its  argument  matrix. 

O  j. 

To  derive  the  filter  as  the  sequential  processor  implemented  in  MINIVAR,  H 
a  regression  formula  will  be  applied  to  incorporate  into  the  estimate  each  new 
data  point  as  it  arrives: 

n/n)  =  y(n/n-l)  +  9y(n) 


z( 


£Tz(n)  -  A^Cn/n-lTj 


(4.8) 
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where 


B^(n)  is  a  9  x  (dim  z)  gain  matrix  which  has  yet  to  be  determined,' 


Az(n/n-l)  =  h  [x(n/n-l)]  -  h  [x°(n)]  in  the  notation  of  (4«3  )• 


The  computation  of  the  one-step  extrapolation  y(n/n-l)  is  straightforward, 

given  state  and  parameter  initial  conditions  (or  prior  in-step  estimates) 

y( n-l/n-1).  From  equations  (4.4),  since  w,  (t)  has  an  expected  value  of  zero, 

a  a  -(t  -t  ,)/T 

A//,\  A  /  ,  /  ,  \  n  n-1  '  a 

u^(n/n-l)  =  u^( n-l/n-1)  e  , 

u2(n/n-l)  =  u^  (n-l/n-1) 

u^(n/n-l)  =  u^  (n-l/n-1)  (4.9) 

The  rest  of  the  y(n/n-l)  vector  is  computed  directly  from  (4.1)  by 
numerically  integrating  up  to  time  t  =  t  after  the  appropriate  rotational 
changes  have  been  made: 


d(Ax)/dt  =  f(x°  +  Ax,  t)  -  f(x°,  4)  +  F(x°  +  Ax,  u,  t) 

Ax  |  =  Ax  (n  -1/n-l),  (4.10) 

initial 
cond . 

where  Ax  and  u  are  short-form  notation  for  Ax(t/t^  1)  and  u(t/t^  ^),  respectively. 

Hence.  (n)  in  (4.8)  is  the  only  computation  for  which  the  machinery  does 
not  yet  exist.  The  remainder  of  the  section  will  describe  how  it  is  computed 
by  linearization  techniques.  Note,  however,  that  up  to  this  point  in  the 
discussion  no  linearizing  assumptions  have  been  made,  since  so  far  B  (n) 

y 

could  just  as  easily  have  been  a  function  of  jr  as  not. 

"';jThe  notation  dim  z  is  the  dimension  of  the  vector  z. 
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B.  MINIMUM-VARIANCE  DERIVATION 


1.  Basic  Principles 

By  their  very  mode  of  operation,  ground-based  sensors  deal  with 
sampled-data  information  which  can  become  rather  sparse  after  normal 
pre-processing  and  conditioning  is  performed.  Hence,  only  at  discrete 
points  in  time  need  we  have  the  ideal  curves  (or  their  equations)  to 
which  the  data  is  to  be  fit. 

This  section  will  generate  the  sampled-data  expressions  for  y  (t) 
which  describe  the  ideal  orbital  motion.  They  will  not  be  integrated  to 
yield  the  best  estimate  itself,  since  equation  (4.10)  has  already  been 
developed  to  do  that.  Rather  these  sampled-data  expressions  will  be 
used  for  the  received  data. 

Because  of  the  inherent  theoretical  problems  in  obtaining  B(n), 

the  equations  of  motion  must  be  linearized  around  a  nominal  trajectory, 

and  Gaussian  probability  distributions  will  be  assumed  for  all  the  random 

6  7  11 

variables.  Then  some  rather  well-developed  theory  ’  ’  and  available 
computer  programs  ^  ^  apply  after  some  modification. 

2.  The  Linearized  Equations  of  Motion 
Consider  equations  (4.1)  and  (4.10): 

d(Ax)/dt  =  f(x,  t)-f(x°,t)  t-F(x,u,t), 
d(Ax)/dt  =  f  (£,t)-f  (x°,  t)  HF(x,u,t) 

where  again  Ax  and  u  are  A*(t/tn  ^  )  and  u(t/tn  ^)  respectively,  and  x(t) 
x°(t)+Ax(t)  and  x(t/tn  . )  =  x°  (t)  +  Ax(t/tn  ^). 

Note  that  if  x  and  $  are  both  sufficiently  close  to  x°  and  hence 
to  each  other,  then  the  f  terms  can  be  replaced  with  the  first-order 
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(4.11) 


o 

terms  in  the  Taylor's  series  expansions  around  x  : 

a  f 

d  ( Ax) / dt  ~  (  ^ Q  Ax  F(x^Ujt)j 
af 

d  (Ax)  /dt  (_^_)o  Ax  F(x,u,t).  (4.12) 

Now  F(x,u,t)  can  be  expanded  in  a  Taylor's  series  around  the  best  estimates 
x  (t/t  ,)  and  u  (t/t  ,):  to  first  order, 

aV  aV 

F  (x,u,t)**  F  (x,u,t)  (■  a  —  )  ( Ax-4x)  (yy — )  (u-u)  (4.13) 

Pf  a  f 

where  ( — - - )  is  ■ —  evaluated  at  x(t/t  ),u(t/t  ). 

v  ax  dx  -v  '  n  —  n' 

Moreover,  we  can  define  the  error  vector 


£a  (t/t^)  =  a  (t)  -  a  (t/t,  )  for  any  vector  a  (t).  In  the  special 
cases  where  t  assumes  discrete  values,  t  ,  the  notation 
will  be  simplified  to  e  (n/k)# 

3. 

Clearly,  then,  equation  (4.12)  can  be  subtracted  from  (4.11),  after  (4.13) 
is  substituted  in,  to  yield  the  error  propagation  equation 


e 

-x 


*  a  f  Pf 

( - )  +  ( — 

v  <3  x  '  o  vdx 


e  (t/t  ,) 
-x  v  1  n-1' 


/■V 

a  F 
a  u 


^  “U  (t/tn_1) 


(4.14) 


A 

d  f  d  F 

In  the  coefficient  matrix  (— - )  i-  (— - )  the  non-Keplenian 

v  d  x  o  dx  * 

accelerations  appear  in  the  lower  three  rows.  However,  the  components 
of  F  are  smaller  than  those  in  f  by  at  least  two  orders  of  magnitude 
even  for  the  lowest-altitude  satellites.^  Since  the  components  of  the 
lower  half  of  f  attenuate  as  l/r  and  those  of  F  attenuate  at  least  as 
fast  in  the  altitude  range  of  interest,  it  is  clear  (see  the  slopes  of 
the  curves  in  Figure  3)  that  the  column  inequalities 
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|  8F±/dx\«  Idf./dXjl  ,  j  =  1,2,3,  all  i, 
hold  in  a  component-by  component  sense  when  averaged  over  a  complete 
orbit.  Because  F  contains  velocity  terms  (xy  Xy  x ^)  whereas  f  does 
not,  the  column  inequalities  will  fail  to  hold  for  j  =  4,  6.  Still, 

the  velocity  contributions  (damping  terms)  are  very  small,  and  a  part 
of  these  are  included  in  (4.14)  via  the  term 


A 

dF 

du 


e 

— u 


Since,  in  addition,  x°(t)  is  forced  by  the  rectification  process 

described  in  Section  III  to  remain  close  to  the  best  approximation  to 

df 

the  actual,  decaying  orbit,  we  can  expect  ( — t - )  to  provide  suf- 

°i  °  aF 

ficiently  accurate  long-term  behavior  in  (4.14)  to  allow  (— — )  to 
be  dropped. 


FIGURE  3*  Relative  Magnitudes  of  Keplerian 
versus  Pertubative  Accelerations 
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All  such  terms ,  except  one,  will  in  fact  be  dropped.  The  one  that  will 
be  kept  is  the  oblateness  contribution.  This  is  because  in 


F(x,  u,  t)  = 

0 

+ 

0 

+ 

0 

+ 

’  0 

/-drag 

-^-solar . 

•  • 

r 

^sun,  moon  _ 

F 

L  — 0  J 

the  oblateness-acceleration  3-vector  is  a  function  of  the  geopotential 
coefficients  J,  which  correspond  to  all  but  the  first  three  entries  in  the 
vector  that  appears  in  (4.14).  However,  the  satellite  does  not  sense 
the  J  vector;  rather  it  senses  the  oblateness  force,  which  happens  to  be 
modeled  often  as  an  expansion  (3*8)  having  coefficients  J.  There  are  only 
three  components  of  force,  but  (in  our  case)  113  components  in  J.  Hence, 
to  account  for  these  modeled  but  not  actively  estimated  parameters,  it  is 
economical  to  let  the  3-vector 


X2(t)  =  F^x,  £)  -  FqCx,  j) 


denote  the  change  in  the  force  due  to  deviations  in  position  and  in  J  from 
their  best  estimates.  Again  employing  the  first-order  terms  in  a  Taylor’s 
series  expansion  about  x,  J,  we  find 

A  A 


v2(t) 


~  (d?)  -  ^(t/tn-l)] 

•A 

dF 

=  (5?)  [^(t)  -  ArCt/t^)] 


dF 

+  >  W, 

A 

6F 

+  <aT>  «. 


(4.15) 


where  r  and  Ar  are  the  first  three  components  (the  position  components)  of 
x  and  Ax,  respectively.  Obviously,  since  v  (t/t^  n )  =  0,  then 

-V  ^An-l)  (4.16) 
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Identifying  the  first  three  components  of  e  in  (4.14)  as  e  and 


-u 


ft 


1 


the  rest  as  e^  ,  we  find  that  after  throwing  away  all  of  (— )  except  the 
oblateness  part  and  using  (4.4), 


e  (t/t  . ) 
-x  '  n-1 


e  (t/t  ,) 
-v^  '  n-1  J 


A 

af 

<ei>o 

aF 

!  ^ 

1  “1 
! 

_ l 

prt' 

0 

0 

0 

e  ( t/t  , ) 
-x  '  n-1 


e  (t/t  , ) 
-v^  '  n-1 


+ 


w-L(t) 


0 

0 


(4.17) 


Even  though  e^  is  a  function  of  e^,  the  equation  can  still  be  integrated 
formally  in  the  following  way.  Write  only  the  first  part  of  (4.17): 


e  (t/t  ,)  =  A(t)  e  (t/t  ,), 
'  n-1  ^y  '  n-1  5 


(4.18) 


where 


e 

-y 


e 

-x 


e 

— v. 


1  J 


A(t)  = 


af 

^ax^o 

!  aF 

W 

_ 

j-1^ 

0 

!  0 

i 

'  0 
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Find  the  transition  matrix  $  (t,  s)  for  (4.18).  This  can  be  done  either 

y 

by  noting  that 


dy(t) 

^y(t *  ~  dy(s) 


dx(t) 

dx(t) 

dx(s) 

H 

dv^s) 

T  7  1  1 

dv1(t) 

dv1(t) 

dx(s) 

dV-^s) 

(4.19) 


or  by  writing  the  Voiterra  series  solution  for  (4.18) , 

t  t  CL 

$y(t,  s)  =  I  +  J  A(a)da  +  f  J  A( a)A(a^)da^da  + 

s  *s  s 


t  a  cl, 


nr  A( a)A(  )A( a2)da2da1da 


s  s  s 


and  picking  out  the  developing  infinite  series. 


If  we  take  the  first  approach  to  finding  J  ,  we  write 

y 


x(t)  =  x[v,(t),  V  (t),  x(s),  t,  s]  , 


-(t-s)/rd 

- 

e 

1 

vx(t)  = 

1 

vx(s)  + 

0 

1 

0 

t 

/ 


-(t-a)/r. 


w^(a)da, 
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and  we  abreviate  equation  (3.4)  to 


x(t)  =  f(x,t)  +  F(x,  u,  t)  =  W  [x(t),  v1( t ) ,  v2(t),  t  ]  .  (4.20) 

Applying  the  chain  rule  to  (4.20) ,  we  obtain 


dx(t)  dW(t)  dx(t) 

dx(s)  dx(t)  dx(s) 


(4.21a) 


dx(t)  dW(t)  dx(t) 

d^TiT  dx(t)  dv^(s) 


dW(t) 

^TiT 


(4.21b) 


dv,(t) 

aitrr"  °> 


(4.21c) 


dv1(t) 

dv^(s) 


"  -(t-s)/rd 

e 


(4.2ld) 


where  W(t)  =  W  £x(t),  v^(t),  v  (t),  .  Again,  to  first-order  accuracy, 


dW 

dx 


df 

<S>o 


A 

dF 


df 


+ 


~  (ta}o 


Hence  (4.2la)  is  a  homogeneous  equation  in 

df 

two-body  coefficient  matrix  (— )q,  so  that 


dx(t) 

dx(  s ) 


with  approximately  the 


dx(t) 
dx(  s ) 


5=3 


$(t,  s). 


(4.22) 
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On  the  other  hand,  (4. 21b)  becomes 


dx(t)  df  dx(  t)  dW(t)  dv^t) 

av,^( s )  ^ax^o  av1  ( s )  4  ^av-^tr  ^  ^(s)) 


ax(t) 

Integrating  —  with  respect  to  time  t,  we  have 


ax(t)  ax(s)  p  _  _ 

aPT?y  ~  $(t’  5)  aVTTJ  H  J  *(t’  u)  (aPTi7)  (o^ti7)  du' 


aw(u)  av1(u) 


-i'  '  -v 


(4.23) 


Since  by  definition 


$y(s,  s)  5  I, 


3x(  s ) 


the  initial  condition  on  (4-23)  is  - — 7 — v  =  0  for  all  s.  Moreover,  from 

avx(s)  5 

(3.4)  and  (3-5),  we  obtain  the  6x3  matrix 


aw(t) 

aiJTtT 


0 

3(t) 


0 


£(t) 


0 


Y(t) 


which  in  conjunction  with  (4-2ld),  provides  the  final  form  for  (4.23): 


ax(t) 

av1(s) 


ax 


z 

(t,  s)«y  $(t,  a) 


§.(a)< 


-(a-s)/r. 


0 


2.(a) 


0 


r(a) 


da, 
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(4.24) 


A  A  A 

where  the  notation  §_(a)  and  Y(a)  indicates  evaluation  around  x(a/s). 


Equations  (4.2lc),  (4.2ld),  (4.22),  and  (4.24)  define  the  transition 


matrix  fy(t,  s)  in  (4.19): 


$y(n,  n-l)  = 


$(n,  n-l) 


0 


dj,  X 


^1 
"-Tn/^d" 


where  T  =  t  -t  n 
n  n  n-l 


The  formal  solution  to  (4-17)  can  now  be 
written  as  the  Volterra  equation 


e 

-y 


(t/t  .)«  I  (t,  t  _  )e  (t  ,/t  .)  +  /  §  (t,  s) 
'  n-l  *y  *  n-l  Y  n-l'  n-l  J  *y  * 


n-l 


/  v*. 


s) 


"n-l 


0 


*‘Ci7 

0 

-0  J 


ds 


(4.25) 


which  in  general  cannot  be  solved  in  closed  form  because  the  middle  term  is 

a  function  of  e  .  However,  that  term  is  extremely  small  compared  to  the 

rest  of  the  expression.  Hence  it  is  not  unreasonable  to  approximate  that 

term  by  rectangular  integration:  that  is,  to  take  e  (s/t  .. )  to  be  constant 

"V2  n 

between  time  points  t  and  t  where  these  instants  may  be  data  points,  or 
rectification  points  called  by  the  Encke  integrator  (see  Section  III),  or 
rectification  points  called  by  a  test  on  the  constancy  of  e^  (s/t^  ^). 
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The  approach  will  be  to  take  e  (s/t  , )  =  e  (t  ,/t  ,)  for 

^  -v  '  '  n-1  -V  n-1'  n-1 

t  t£t  Proc*uce  staircase  approximation  suggested  in  one  dimension 
by  Figure  4,  below.  Of  course,  we  cannot  test  the  error  itself  to  ascertain 
that  it  is  near  enough  to  being  constant,  but  we  can  test  to  see  if  its 
covariance  has  changed  significantly  with  reference  to  the  covariance  of  that 
to  which  it  contributes:  i.e.,  the  ®y( tAn_i )  in  equation  (4.25).  This  will 
be  made  clear  shortly. 


FIGURE  4.  STAIRCASE  APPROXIMATION 

TO  e  (t/t  -  ) . 

-v  '  n-1 


A  final  linearization  is  still  to  be  performed  on  the  observation 
equations.  In  particular,  using  the  notation  developed  in  this  perfect-sensor 
case  for  (4.2)  and  (4.8),  we  see  that  the  residual  term  in  (4.8)  upon  which 
B(n)  operates  is 


Az(n)  -  Az(n/n-l)  =  hjjc(n)j  -  -  h  [x(n/n-lj]  +  h^sAn^  , 


which,  after  a  first-order  expansion  around  x(n/n-l),  can  be  written  as 

A 

dh 


Az(n)  -  Az(n/n-l)*s(— )  [x(n)  -  x(n/n-l)J 

=  H(n)  [Ax(n)  ~  Ax(n/n-l)J  , 


(4.26) 
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dll 


where  H(n)  =  (— )  is  the  partial  derivative  evaluated  at  x(n/n-l) 


3 .  Prediction  Covariance  Computations 

Evaluate  equation  (4*25)  at  time  t  =  t 

Call 


t 

1  rn  , 

—  (n,  n-l)  =  /  <|  (t  ,  s) 

I2  J  y  n 


6v 


ds 


n-l 


0  0  0* 
0  0  0 
0  0  0 
10  0 
0  1 


0 

0  0  1 
0  0  0 
0  0  0 
0  0  0 


(4.27) 


and  assume  the  staircase  approximation 


e  (s/t  )=e  ( n-l/ n-l). 

-v  '  n-l  — ' v  '  ' 

2  2 


(4.28) 


Define  the  9-vector 


t, 

w(n-l)  =J  Ij/V  s) 


"n-l 


w1(s) 

0 

0 


ds , 


■/' 


"n-l 


dx 

du^  ^  tn*  ^ 

0 

0 


w1(s)  ds, 


(4.29) 
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and  its  covariance  matrix 


Cov  cj(n)  = 


Vn) 


«u(n> 

!  q(n)  0  I  0 

'  —  •  * 

qT(n) 

I  «  (n) 

■  -  1 

0 

i  0 

-  -i 

0 

1 

1  0 

where,  from  (4.29)  and  (4.5),  it  can  be  computed  that 


Qn(n-1) 


q(n-l) 


<r2(n-l) 


-  < 2/  V  a  drag 


f  &£(tn*s)  {s£(Vs)}Td= 


t  , 
n-1 

r  t 


fn  4*<t 

)  du,  n 

t  1 

tn-l 


XT*  (t  ,s)e 


-(t  -s)/r. 


ds 


-2(t  -t  ,  )/r 
-(l-e  n  n-l/Cd)(j2 


drag 


Then  (4.25)  can  be  rewritten  as 


aZ 

e^n/n-l)  =  $y(n,  n-l)  e^n-l/n-l)  +  —  (n,  n-l)  ey  (n-l/n-l) 


+  to  ( n-l ) , 


(4.30) 


(4-31) 
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and  its  covariance  P  (n/n-l)  becomes,  after  we  note  that  w(n-l)  is 

yy 

uncorrelated  with  co(n-2),  and  hence  is  uncorrelated  with  everything  else 
on  the  right, 


Vn/n-1)  =  Jy(n,  n-1 ) Pyy( n-1/ n-1  )$y4 „ ,  n-1) 

+  ®y(n,  n-ljp^^tn-l/n-l)  (n,  n-1) 

+  f|-  (n,  n-DP^  T(n-l/n-l)®  T(n,  n-1) 

2  J 

^  T 

+  (n,  n-l)Pv^(n-l/n-l)  n_i)  +  (^(n-l)  (4.32) 


A  recursion  can  be  developed  once  the  in-step  covariance  matrices 
P  (n/n),  P  (n/n),  and  P  (n/n)  are  defined  in  terms  of  covariance 

yy 7  2  V2V2 

matrices  with  argument  (n/n-l).  However,  since  these  in-step  covariances  are 
functions  of  B  (n),  that  gain  matrix  will  now  be  developed  explicitly. 

y 


4.  Computation  of  Optimal  Gain  B  (n) 

y 

We  begin  by  applying  the  linearization  (4-26)  to  (4.8)  and  subtracting 
y(n)  from  both  sides  of  the  resultant  equation.  We  find 


e  (n/n)  =  e  (n/n-l)  -  B  (n)H(n)e  (n/n-l) .  (4-33) 

u  v  « y  ^ 

If  we  define  the  (dim  z)  x  9  observation  matrix 

Hy(n)  =  [H(n)  o],  (4-34) 
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then  (4.33)  becomes 


®y( n/n )  =  (I  -  ByHy)ey(n/n-l)J  (4-35) 

where  the  argument  (n)  has  been  dropped  from  B  and  H  . 

y  y 

Clearly,  then, 

P  (n/n)  =  ( I-B  H  )  P  (n/n-l)  (I-B  H  f ,  (4-36) 

yy  y  y  yy  y  y ’ 

and  its  total  differential,  given  that  B  is  the  only  variable  that  can 

y 

be  manipulated,  is  just 

dP  (n/n)  =  dB  I[h  P  (n/n-l )H  T~|  B  -HP  (n/n-l )\ 

yy  y  U  y  yy  y  J  y  y  yy  J 

+  dB  f  [h  P  (n/n-l)H  T]  B  ‘  -  H  P  (n/n-l)}  x  (4-37) 

y  1 1  y  yy  y  1  y  y  yy  J 

Now,  to  achieve  the  minimum  required  by  (4*7)?  the  total  differential  of 
the  terms  on  the  left  of  that  expression  must  be  zero: 

d  [trace  {AP  (n/n)A)^  =  trace  { AdP  (n/n)A)  =  0.  (4-38) 

yy  yy 

It  is  a  property  of  the  trace  that 

trace  {xy}  =  trace  $YX)  , 

traceJx+Y^  =  trace  X  +  trace  Y, 

T 

trace  X  =  trace  X  , 

for  any  square  matrices  X  and  Y.  Therefore  (4-37)  and  (4-38)  combine  to 
yield 

trace{A2{By[HyPw(n/n-l)HyT]  -  n/n-l )HyT )  dBy}  -  0. 
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Since  this  must  hold  for  all  differentials  dB  , 

7 


A2(b  [h  P  (n/n-l)H  T  1  -  P  (n/n-l)H  T  )  =  0  (4-39) 

'  7  1  7  yy  7  J  yy  7  ' 


is  the  condition  for  an  optimum.  When  A  is  any  nonsingular  matrix , 
the  optimal  B  (n)  is  defined  to  be 

y 


B  (n)  =  P  (n/n-1)  H  T(n)  [H  <n)P  <n/n-l)H  T(n)] 


yy 


yy 


(4.40) 


6  7 

which  is  the  usual  Kalman  gain  matrix.  5' 

In  the  event  that  we  are  interested  only  in  estimating  the  mass  parameters , 
we  can  set  all  the  cu  in  A  to  zero  excep  for  ag  and  .  Then  only  the  last 
two  rows  in  B  (n)  are  defined  by  (4*39),  the  others  being  arbitrary.  These 

y 

last  two  rows  are  the  same  as  the  last  two  rows  in  the  usual  Kalman  matrix 
given  by  (4.40).  The  remaining  rows  we  can  set  to  anything.  Let  us  set  them 
to  the  usual  Kalman  values  defined  by  (4.4Q).  Hence,  independent  of  the 
weights  A,  the  Kalman  gain  (4*40)  is  always  optimal .  Note  that  if  we  had 
generalized  A  to  be  non-diagonal,  the  result  would  have  been  the  same. 


In  terms  of  the  partition 


P  (n/k)  = 

yy 


!  Pxv1(l'/k) 

P^’tnA)  |  Pv  ("A)_ 


equation  (4.40)  becomes 


By(n)  = 


\(n) 

P  (n/n-1) 

XX 

A(n). 

P  1 (n/n-1 ) 

L^i  J 

HT(n)  [H(n)Pxx(n/n-l)HT(n)]  (4.4o’) 
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where  B  (n)  is  the  6  x  (dim  z)  gain  for  the  satellite  state  estimates  and 

B  (n)  is  the  3  x  (dim  z)  gain  for  the  mass  parameters  v.. . 

V1  1 


5 .  In-Step  Covariance  Computations 

Inserting  (4. 4-0)  in  (4.36),  we  can  show  that 


P  (n/n)  =  ( I-B  H  )2P  (n/n-l)  =  ( I-B  H  )  P  (n/n-l).  (4.41) 

yy  yy  yy  yy  yy 

T 

Moreover,  we  can  post-multiply  equation  (4.35)  by  (n/n),  use  (4.15) 
and  (4.16)^  and  take  the  expectation: 

P  (n/n)  =  (I-B  H  )  P  (n/n-l) 

yv.  y  y  yy 

+  (I-B  H  )  P  .( n/n-l) 
y  y  yJ 

where  the  argument  (n)  on  the  partial  derivatives  indicate  evaluation  at  the 
point  x( n/n ) ,  and 


multiplying  (4-31)  by  (aJ)T  and  taking  the 

expectation: 


dF 

-o 

dy 


dF 

— o 

dr 


000  | 000 
000  1 000 
000  1000 


The  matrix  P  .( n/n-l)  is  found  by 


A  <T> 

dF 

5T2  (n) 
dy 


A  <71 

dF 

— o 

dJ 


(n), 


(4.42) 
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r  6F 

Pyj(n/n-l)  -  U(n,  n-1)  +  —  (n,  n-1)  (n-1) 


PyJ(n-l/n-l) 


A 

dI  dF 

+  3E'n>  »-1)af(n-1)  PJJ- 


where  Pjj  =  Cov  (aJ).  The  in-step  matrix  T(n-l/n-l)  obtains  by 


yJ 


multiplying  (4*35)  by  (AJ)X  and  taking  the  expectation: 


P  T(n/n)  =  ( I-B  H  )  P  T(n/n-l). 

yJ  '  y  y  yJ  ' 


This  completes  the  detail  for  (4.42). 
so  far, 

^(n/n)  =  D(n)P^r(n/n-l), 


To  summarize  what  we  have  done 


PyT(n/n)  =  D(n)  PyJ(n/n-l), 

T  T 

dF  6F 

P  (n/n)  =  P  (n/n)  (n)  +  P  T(n/n)  ~rr—  (n), 

y^2  yy  sy  yJ  aJ 


where 


D(n)  = 


”l-Bx(n)H(n)  | 

i 

0 

l 

-B  (n)H(n)  | 

I 

and  where  P  T(n/n-l)  is  given  by  (4-43) 
y<J 


(4.43) 

(4.44) 

(4.45) 

(4.46) 

(4.47) 

(4*48) 
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To  compute  P  (n/n)  from  (4.15)  is  now  theoretically  straightforward: 
V2V2 

A  A«t* 

dF  dF 

P  (n/n)  =  ~~  (n)P  (n/n)  (n) 

v  Vg  dr  rrv  '  '  dr  v 


.a  /\  T 

dF  t 

+  tt-  (n)Pjj  —  (n)  +n(n/n)  +  TT  (n/n), 


dJ 


(4.49) 


where 


TT(n/n) 


A  At* 

dF  dF 

S?  (n)  PrJ(n/n)  W  (n)’ 


and  P  and  P  T  are  the  submatrices  defined  by 
rr  rJ 


P  (n/n) 
xx  ' 


P  (n/n) 
rr  ' 


P  -(n/n) 
rr 


P.  .(n/n) 
rr  ' 


PyJ(n/n) 


prJ(n/n) 


P.j(n/n) 


pViJ(t./t.) 

the  latter  being  given  by  (4*46). 

3o  far,  the  updating  of  P  has  been  described  only  for  the  cas'e 

yy 

of  closely  spaced  data  points,  t  .  and  t^;  for  widely  spaced  points, 
the  assumption  of  constant  (sAn  q_)  in  equation  (4.25)  is  not  valid. 
This  problem  is  easily  circumvented;  there  is  no  need  to  wait  until  data 


point  t  to  perform  the  covariance  matrix  update.  Rather,  when  (accord¬ 
ing  to  some  criterion)  the  assumption  of  constant  e  is  in  danger  of 

V2 
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being  violated,  say  at  intermediate  time  t  ,  ,  =  t  ,  +  At,  a  new  in- 

n- 1  ,  1  n- 1 

step  covariance  matrix  can  be  defined  as 

Pyy^tn-l,l/tn-l,l)  =  Pyy(tn-l,l/n)> 

i.e.,  the  filter  gain  B  (t  ,  , )  that  appears  in  equation  (4.41) 

7  n--L,  -L 

is  set  to  zero.  The  update  procedure  for  the  extrapolated  covariance 

matrix  can  now  proceed  normally  from  tn_^  either  to  the  next  data  point, 

t^,  or  another  non-data  update  point,  tR  ^  2 • 

Calculate  P  (n-l/n-l)  from  (4.49);  compute  P  (t  ,  ,/n-l)  from 
v2v2v  '  y  yy  n-l,r 

(4.32),  retaining  the  sum  of  the  first  and  last  terms,  which  involve  P  (n-l/n-l) 

yy 

and  Q  (n-l).  Compute  now,  the  matrices  P  (t  ..  ,/t  ,  ,  )  and 

yy  F  ’  yv2v  n-l,l'  n-l,i' 

P  (t  ,  ,/t  ,  ,  )  to  be  used  in  the  next  update.  Before  updating, 

n- l ,x  n— x,x 

recompute  Pyy(’*:'n_p  ^/n-l)  using  the  same  first  and  last  terms  as  before. 

In  place  of  P  (n-l/n-l)  and  P  (n-l/n-l),  use  P  (t  ,  ,/t  ,  , ) 

y  yv2v  '  '  v2v2v  '  '  yv2  n-l,  r  n-l,  1 

and  P  (t  ,  ,/t  ,  ,),  respectively. 

n-1,1  n-l, 1  y  * 

The  changes  in  the  P  and  P  matrices  are  indicative  of  the  errors 

yv  v„v„ 

J  2  2  2 

implicit  in  the  constant  e  assumption;  the  difference  between  the 

V2 

P  (t  ,  ../n-l)  matrices  using  the  old  and  new  P  and  P  illustrate 
yyv  n-i,r  &  yv2 

the  order  of  magnitude  of  the  maximum  errors  in  P  that  result  from  the 

yy 

assumption  under  test.  Thus,  for  the  elements  p. .  of  the  extrapolated 
1 

P  ,  and  p. .  of  the  recomputed  P  ,  find 

yyJ  *  yy 

>  ib.j--.Xi 


Then,  if  E  >b,  for  some  upper  bound,  b,  reduce  the  At  used  to 

yy 

advance  to  the  next  non-data  update  point,  if  any.  If  E  <a,  for  some 

yy 

lower  bound,  a,  increase  At. 
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C.  SUMMARY  OF  FILTERING  PROCEDURE 


1.  Start  with  an  initial  estimate 


z  (o/o) 


Ax(0/0) 


l2]_(o/o) 


of  the  satellite  states  &x  and  the  mass  parameters  v^. 

Use  the  mass-parameter  model  (A. 9)  to  compute  the  extrapolations 


Zi  (Vo)  = 


~-(t-t  )/r 
e  o"  d 


v  (0/0). 


3 .  Employ  these  extrapolations  in  the  Encke  equation 

[equation  (4.10)] 

d(Ax)dt  =  f(x°+Ax,t)-f(x°,t) 

+  F(x°+£a,  v^t/O),  J), 

and  integrate  from  the  initial  condition  Ax(0/0)  to  obtain 
A^(l/0).  Also  evaluate  (t/0)  at  t  =  t.  to  obtain  v-,(l/0) 

4.  Guess  at  the  covariances  on  the  initial  estimates: 

V°/°)>  Pyv  <°/°>-  PyJ(°/0)- 

5.  Use  these  matrices  to  compute  [from  (4.49)] 

ss  /v  m 

dF  dF 

P  (0/0)  =  ~  (0)P  (0/0)  (0)  + 

V^v0v  '  '  dr  '  '  rr  '  '  dr  '  ' 

A  AT 

dF  er  t 

+  ^  (0)  Pu  ^  (0)  +  a  (0/0)  +  E  (0/0), 


JJ  dJ 
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where 


/N 

a4 


A  m 

«F. 


a (0/0)  «^r  (0)P„.,(0/0)  (0). 


rJv 


dJ 


6.  Then  evaluate  (4.32): 


Pyy(l/0)  =  ^(1,0)  Pyy(0/0)  |yT(l,0)  +  A(  l/O )  + 

T 

+  at(i/o)  +  d_l_  (1,0)  P  v  (0/0)  (] 


^  V  V 

dI2  2  2  dv2 


V(0) 


where 


A(i/o)  =  $J1,0)  p„r  (0/0)  —  (1,0). 


yv- 


Also  evalute  (4-43): 


pyJ(i/o)  - 


djr  0F 

|  (1,0)  +  v —  (1,0)  —  (0) 

*y  dZ 


pyJ(0/0) 


<3^;  OF 

+  av;  <M)  a?  <°>p 


JJ. 


7-  Then  use  the  first  of  these  in  (4-40  ): 


B  (1)  = 
y 


"b  (1) ' 

X 

[V1/0)' 

t  r 

B  (1) 

V1  J 

-Pxv<1/°2 

HT(1) 

IKD 


8.  Use  B  (1)  in  (4-45): 

<J 


Pyyd/l)  -D(l)  Pyyd/O), 


,0)  + 


Pxx(l/0)H  (1) 
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and  (4.4 6): 


VlA>  ■  B(i)pyJ(1/0)> 


where 


D(l) 


I-Bx(l)H(l)  j  0 


-B  (l)H(l)  !  I 

1  i 


9.  From  Pyy(Vl)  compute  P  (l/l)  using  (4.47), 

SS  rn  /\  m 

6F  dF 

P  (1/1)  =  P  (l/l)  ~~ —  (1)  1  D(1)P  r(l/0)  (1), 

yv2v  '  yy  ay  yJ  aj 

and  compute  P  (l/l)  from  step  5  with  "0"s  replaced  by  "l"s. 

V2V2 

10.  Use  B^(l)  and  y(l/0),  together  with  a  measurement  z.(l),  to 
compute  y(l/l)  from  (4.8): 

y(  1/1 )  =  y(l/0)  t  B(l)  [  Z(l)-h  [x(l/0)]] 

11,  Return  to  step  1  and  re-do  the  process  with  n0M  subscripts  and 

arguments  replaced  with  nlMs  and  ,!lMs  replaced  with  M2ns.  Skip 

steps  k  and  5.  Repeat  for  t^?  t^,  .  ..,tn, 

The  above  holds  for  two  data  times  t  and  t  . n  sufficiently 

n  n+1  J 

close.  If  the  test  on  the  constancy  of  e  (t/tn)  requires  an  update 

V2 

point  at  tn+At  <  t  the  same  process  as  the  above  is  followed, 
except  that 

(a)  tn+At  replaces  t^^  in  the  program; 

(b)  step  7  is  replaced  with  a  step  that  sets  B  (l)  =  0. 

y 
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D.  SENSITIVITY  COEFFICIENTS 


One  of  the  major  aims  of  this  study  is  to  determine  the  quantitative 
effects  of  an- error  source  on  the  estimates  of  the  ballistic  coefficients 
u^  and  These  effects  are  computed  classically  by  assuming  that  the 

error  contributions  from  each  source  are  small.  In  addition,  it  is  assumed 
that  the  functional  relationships  between  the  error-source  parameters, 
say  a  =  column  (  a.  ,  a_  )  and  the  ballistic-coefficient  estimates 
are  known  to  be,  say, 

£3  =  8j(  E  )  •  (4.50) 

Then  a  Taylor's  series  expansion  in  the  errors  can  be  performed  and 
truncated  at  the  first-order  terms  to  yield 

T 

A  u2  =  (grad  g2)  A  a  , 

A  u^  =  (grad  g^)T  A  a  ,  (4-51) 

where  the  gradients  are  taken  with  respect  to  and  evaluated  around  the 
nominal  values  of  the  parameters.  The  biases  in  the  ballistic  coefficients 
are 

E(A  u2)  =  (grad  g2)T  E(a  a  ), 

E(A  u^)  =  (grad  g^)1  E(A  a  ). 

The  variances  are 

(grad  g2)T  E(A  a  A  a  T)(grad  g2)  -  E^(a  u2), 

(grad^)T  E(a  a  a  a  T)(grad  g^)  -  E2(a  u^), 
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or 


'  =  (grad  g2)  Cov  (a  a  )  (grad 


rp 

=  (grad  g^)  Cov  (a  a  )  (grad  g^). 


(4.52) 


In  terms  of  these  second-order  statistics,  the  gradient  vectors  define 
completely  the  importance  of  each  error-source  parameter  in  the  ballistic- 
coefficient  error  build-up.  The  elements  of  these  gradient  vectors  are 
often  called  "sensitivity  coefficients." 

Unfortunately,  in  a  maximum-likelihood  or  minimum- variance  estimation 
procedure  of  the  sort  required  here,  the  functions  g.  and  gq  in  (4-50) 
cannot  be  written  as  closed-form  expressions  of  the  error-source  parameters 
a  .  Hence  the  gradients  cannot  be  obtained  analytically. 

1 .  Modeled  Parameters 


In  the  minimum-variance  estimation  equations  developed  in  parts 


A  and  B  of  this  section,  certain  parameters  have  variance-covariance  values 
associated  with  them.  It  is  said  that  the  errors  in  these  parameters  are 
thereby  "modeled."  The  stochastic  component  of  drag,  u.  ,  and  the  geopo¬ 


tential  perturbation  acceleration  v^  are  the  only  two  modeled  error  sources 


considered  in  this  perfect-sensor  phase  of  the  study.  Call  these  modeled 


error  parameters  & 2’  •••>  anc*  varaances 


Then  the  minimum- variance  equations  yield  implicitly  the  in-step  esti¬ 
mation  errors 
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i  u2  o2  (  <rs^  <rs^  ...)> 

4u3"G3  (  *b2,  ...)• 

To  ascertain  the  importance  of  reducing  the  ignorance  in  one  modeled  para¬ 
meter  relative  to  another,  assuming  small  reduction  in  such  ignorance,  we 
need  to  find  for  all  i 


dtf2  /dor  ~  AO2  /AT 

U2  3i  U2  3i 


d<r2  /dT  »  AO2  /A  Q 
u3  Pi  u2  3i 


where  A  indicates  a  (small)  change  in  the  variable  that  follows  it. 

Each  of  these  partials  is  a  form  of  sensitivity  coefficient  which 
is  really  the  desired  end  product  of  the  study.  The  manner  of  computing 

it  consists  simply  of  making  a  run  with  all  the  reduced  by  the 

*1 

amount  A  . 

pi 

2.  Unmodeled  Parameters 

The  minimum-variance  expressions  contain  a  few  parameters  which 
may  have  incorrect  values,  but  for  which  no  reasonable  variance-covariance 
information  exists.  Such  parameters,  having  so-called  "unmodeled"  errors, 
are  typified  by  the  correlation  time  constant  in  the  stochastic  drag 

expression  (4.4),  or  the  geomagnetic  index  a^  that  enters  the  nominal 


-  45  - 


atmospheric-density  calculations.  Call  these  unmodeled  parameters  ^  , 

&  ,  . . . ,  and  their  (unknown)  variances  cri,  ,  ,  ... 

2  Y1  y2 

Now  the  estimation  errors  will  not  directly  reflect  any  improvement 
in  our  knowledge  of  the  error-source  parameter  values.  However  the  esti¬ 
mates  themselves  are  implicit  functions  of  the  parameters:  from  (4.50), 


U2 

^2  ^1  9  ^  2’  *  *  *  *  ^l9 

' ' *  b 

/s 

u3 

~~  (3 j  ^  2. 9  *  *  * 9  ^1 9 

r  ,  ...). 

Hence  by  changing  each  Y.  while  holding  the  other  parameters  at  their 
nominal  values ,  each  sensitivity  coefficient  in  grad  gr  and  grad  g^  can 
be  determined  numerically.  Put  into  (4.52) ,  these  coefficients  lead  at 
once  to  the  pinpointing  of  the  dominant  error  sources. 

E.  MINIMUM-Y AH IANCE  VERSUS  MAXIMUM  LIKELIHOOD 

Although  the  contract  specifies  analysis  of  the  errors  in  the 
maximum-likelihood  estimation  of  mass,  the  approach  contained  herein  is 
based  upon  minimum- variance  estimation.  As  is  sho\m  in  detail  in  Appendix  V, 
the  theoretical  details  may  be  different,  but  the  practical  solution  to  the 
problem  is  essentially  the  same  in  both  cases.  Computationally,  the 
problems  that  arise  in  one  approach  arise  also  in  the  other. 
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SECTION  V 


ENVIRONMENTAL  MODELS 

A.  SOLAR  RADIATION  PRESSURE 

The  model  is  discussed  in  Appendix  I.  It  will  become  important  at 
altitudes  of  about  500  km  or  higher,  and  will  be  the  sole  mass-determining 
force  at  altitudes  in  excess  of  1000  km. 

B.  ELECTROMAGNETIC  DRAG 

A  detailed  analysis  of  the  three  types  of  significant  electromagnetic 

21 

drag  is  carried  out  in  Appendix  H  based  upon  the  CIRA  1965  tables  for  ion 
and  neutral-particle  densities. 

Of  the  three  types  of  electromagnetic  drag,  two  are  distinguishable 
from  atmospheric  drag  because  they  are  latitude- longitude  dependent  and 
have  a  velocity  dependence  different  from  atmospheric  drag.  These  consist 
of 

(1)  Electromotive  drag,  where  the  charged  vehicle  interacts  with 
the  earth's  magnetic  field; 

(2)  Induced  drag,  where  the  potential  induced  in  the  moving 
satellite  (conductor)  by  the  earth's  magnetic  field  sets  up 
a  current  through  the  ionized  atmosphere  which  produces  a 
back-EMF  retardation. 

Under  the  worst-case  conditions  (the  largest  satellite,  highest  charged- 
particle  densities,  and  stro-ngest  magnetic  field),  these  effects  are  at 
least  two  orders  of  magnitude  smaller  than  atmospheric  drag. 

The  remaining  type  of  electromagnetic  drag  is  indistinguishable 
from  atmospheric  drag  to  the  ground-based  observer.  It  is  known  as  coulomb 
drag,  and  results  because  the  charged  satellite  electrostatically  attracts 
to  it  ionized  particles  that  it  would  not  otherwise  hit.  The  ensuing  increase 
in  the  effective  area-to-mass  ratio  may  reach  10  per  cent  or  more.  It  will  be 
accounted  for  by  making  off-line  corrections  to  the  estimate  of  ballistic 
coefficient 
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C.  GEOPOTENTIAL  MODEL 


The  earth's  geopotential  model  is  discussed  in  Appendix  HI. 

D .  SOLAR-LUNAR  GRAVITATION 

The  solar-lunar  gravitational  accelerations  on  the  satellite  will  be 
computed  from  the  assumption  that  the  sun  and  the  moon  are  point  masses. 

q  l 

Computational  details  will  be  as  specified  for  the  MINIVAR  computer  program,'' 
and  presented  in  equation  (3*7)  of  this  report.  The  respective  distances 
and  gravitational  constants  will  be  assumed  exact . 

E.  ATMOSPHERIC  DRAG 

Atmospheric  drag  can  be  computed  interchangeably  from  the  expression 
of  Appendix  II or  from  the  usual  velocity-squared  law  given  in  (3-5).  Both 
require  an  accurate  upper-atmosphere  model.  The  Jacchia  1965  model,  published 
in  Reference  8  and  described  briefly  below,  will  be  used  here,  in  conjunction 
with  expression  (3*5)  of  this  report.  RMS  values  and  correlation 

times  for  equation  (3-5)  will  be  chosen  in  concert  with  the  contracting 
agency.  Suggested  values  are  given  in  Reference  7- 

The  following  summary  of  the  Jacchia  1965  model  was  extracted,  with 
minor  corrections,  from  Reference  28: 

Jacchia' s  1965  atmospheric  model  (Reference  8^  begins  at  a 
boundary  altitude  of  120  kilometers  where  the  following  assumptions 
are  made : 

1)  =  355°  K 

n(N2 )  =  4.0  x  1011  molecules /cm3 

n(0a )  =  7.5  x  lO10 
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n(0)  =  7.6  x  lO10  molecules /cm3 

n(He)  =  3.4  x  107 

Below  the  altitude,  Z  =  500  kilometers,  n(H)  c  0. 

At  500  kilometers,  the  boundary  condition  for  hydro¬ 
gen  is  : 

LOG10[n(H)]  =  73.13  -  39.40  L0Glo  Ofeco)  +  5  .sClOG^Osoo)]2 


Where ; 

Ten  -  temperature  in  degrees  Kelvin  at  500  kilometers, 

n  =  number  density  of  individual  constituents  of 

the  atmosphere;  nitrogen(N2  )  ,  oxygen(02 ) ,  free 
oxygen(O) ,  helium(He)  and  hydrogen(K) . 

Using  the  boundary  conditions  as  a  starting  point  for 

the  concentrations  n.  of  each  consitituent  i,  the 

1 

following  diffusion  equation  is  integrated  to  find 
the  concentrations  as  a  function  of  altitude  Z. 


dn . 


2) 


n . 

i 


i  =-dZ  _  dT  (1  +  a) 


H. 

i 


Where  : 

T  =  temperature  in  degrees  Kelvin 

Ct  =  thermal  diffusion  factor 

(-  0.38  for  helium 

0.0  for  other  cons i tutuents 
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mass  of  each  constituents. 


m(N3) 

-  4.6515 

X 

o 

1 

grams /molecule 

m(Qa) 

=  5.3129 

X 

lO”33 

grams /molecule 

m(0) 

-  2.6565 

X 

10_33 

grams/atom 

m(He) 

=  0.6648 

X 

lO-33. 

grams/atom 

m(H) 

=  0.1674 

X 

10'33 

grams / atom 

k  =  1.38054  X  10  B  ergs/°K  (Bol tziaann 1  s  constant) 

R  *  6.3  5  6  7  7  X  108  cm.  (radius  of  the  earth) 

2  -o  <  o 

g  =  980.665  (1  +  ~)  cm. /sec.  (acceleration  due  to 

gravity) 


H.  =  -  (scale  height) 

i  m.g 


Equation  2)  when  integrated  becomes: 

Z 


3) 


T^-l  +  a 


n .  = 
i 


nia°  L  T  J 


exp 


"  mi  r  dz"> 

g  T ! 


L  k 


J 

120 


Simpson's  rule  of  order  2  with  a  step  size  of  .1  k.n 
is  used  for  numerical  integration. 

The  following  sequence  of  operations  determine  the 
temperature  T  at  the  desired  altitude  Z 

4)  T  =  -  (T,  -  Tmo)  exp  [-S(Z  -  120)] 

Where : 

T  -  exospheric  temperature  °K 
=  temperature  at  120  km. 
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Z  “  altitude  in  km. 

-X2 

S  -  0.0291  exp  <-y-) 

,  T  -  800.0 
00 

X  =  750.0  +  1.722  X  10J<*  (Tm  -  800. 0)L 

Note:  The  expressions  S  and  X  are  approximations  to  the 
profile  derived  by  L.  Jacchia  by  trial  and  error 
me  thods  . 

5)  P  -T  n^m  gm/cm^ 
i 

p  *  density  as  a  function  of  altitude 

The  following  procedure  is  used  to  determine  the  exospheric 
temperature, 


1)  T  =  357°K  +  3.6  F10.7 

o 

Where : 

T^  is  the  average  nighttime  minimum  temperature 
F10#7  is  the  smooth  solar  flux  over  3  solar  rota¬ 
tions 

2)  ToX  *=  Tq  +  1.8°  (F1c.7  -  Pio.7.) 

Where : 

T^1  is  the  variation  expected  during  a  given 
solar  rotation  (27  days) 

F10.7  is  the  solar  flux  for  the  previous  day. 
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3) 


j^O  .37  +  0.14  sin  (jn  ~ ^)j  F10.v  sin  ^tt 


d  -  5  9^, 
365^ 


Where : 

T  accounts  for  the  semi-annual  variation, 
o 

d  is  the  number  of  days  from  January  1. 


4)  T 


Where : 

T  =  T  (1.0  +  R  cos"1  r\ ) 
D  o 

T  =  T  (1.0  +  R  sinm  0) 
N  o 

r\  =0.5  1 0  -  6  | 

S 

6  =  0.5  |0  +  6  | 

s 


■  5S  ■ 

s  s 

*  ■  C  7FV?) 

V  V 

R  =  0.28 
m  =2.5 

T  =  H  +  (3  +  P  sin  (H  +  a) 

H  TAH_1  (/)-  IAN'1  (^) 
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0  =  -45° 

P  =  12a 

a  =  45° 

T  accounts  for  the  diurnal  effect. 

is  the  daytime  maxima. 

Tvt  is  the  nighttime  minima. 

N 

R  is  a  factor  for  computing  the  maximum  temperature 

as  a  function  of  the  global  minimum  temperature. 

0  is  the  geographic  latitude  of  the  vehicle. 

6  is  the  declination  of  the  sun. 
s 

X  .  Y  ,  Z  are  the  current  coordinates  of  the  vehicle 

V  V  V 

j 

in  an  earth-centered  cartesian  coordinate 
system. 

X  ,  Y  ,  Z  current  coordinates  of  the  sun 
°  v  s 

is  the  hour  angle  of  the  sun 

=  T  +  1.0°  a  +  125°  [l.O  -  exp  (-0.08  a  )] 

P  P 

Where : 

is  the  exospheric  temperature. 

is  the  3  hour  geomagnetic  index  (measured  approximately 
6  hours  before  the  time  in  question.) 


H 


5)  T 

/  rt 
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APPENDIX  I 


PERTURBATIONS  ON  SPHERICAL  SATELLITES  DUE  TO 

SOLAR  RADIATION  PRESSURE 

A.  GENERAL 

A  photon  flux  of  irradiance  I  directed  along  a  unit 

A 

rector  1  from  &  distant  radiating  body  (e.  g.,  the  sun;  to  a 
ss 

satellite  will  produce  a  pressure  of 


e  -  <V°> 


(i.i) 


where  c  is  the  speed  of  light.  When  £  impinges  upon  a  differen¬ 
tial  surface  area  dA,  the  differential  force  due  to  the  incident 
flux  is 


<B  •  <*A>  -  (Vc)  (X-2) 


while  the  differential  force  due  to  the  flux  being  reflected  back 
from  the  surface  is 


(1.3) 


A 

where  i  is  the  unit  rector  in  the  direction  of  the  reflec 
refl 

flux  and  I  is  the  irradiance  in  that  direction, 
r 

The  net  solar  force  due  to  direct  solar  radiation  is  the 
integral  of  (l*j>)  plus  (1.3): 
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where  A  is  the  illuminated  surface  area  of  the  satellite.  For 
s 

satellites  with  uniform  reflectivity,  equation  (1.4)  can  be 
simplified  to  a  form  such  that  the  acceleration  magnitude  due  to  direct 
solar  pressure  is 

direct  "  Fd/m  =  K(Is/c)  (A/m),  (1.5) 

where  K  is  a  parameter  dependent  upon  the  shape  of  the  satellite, 

its  reflectivity  coefficient,  and  its  mode  of  reflection  (diffuse 

on  specular),  A  is  the  cross  sectional  area  as  seen  by  the  radiation 

source,  and  m  is  the  mass  of  the  satellite. 

The  constant  c  is  known  very  accurately,  but  I  ,  which  has 

s 

2 

a  nominal  value  of  1.94  cal/(cm  -  min),*  may  vary  by  as  much  as 

TO 

1.5%  due  to  solar  activity.  ,  'Furthermore,  the  varying  solar 
distance  throughout  the  year  causes  the  irradiation  to  change  by 
as  much  as  3- 5%  from  the  mean.  This  latter  phenomenon  can  be  cor¬ 
rected  for  by  expressing  the  solar  constant  as 

I  -  (R  /r  )2I  (i.6) 

s  v  es'  es'  nom  1  ,D; 

■where  R  is  the  mean  earth-sun  distance,  r  is  the  actual  earth- 
es  es 

o 

sun  distance  and  I  a  is  1.94  cal/ ( cn-min) .  This  leaves  K(A/m), 
which  can  be  considered  a  (partially)  unknown  parameter.  The 
accuracy  to  which  it  can  be  determined  is  highly  dependent  upon 
the  accuracy  to  which  K  can  be  determined  by  optical  means. 

For  non-spherical  satellites,  the  problem  is  considerably 
more  complicated.  This  is  principally  a  result  of  a  different 

*  A  calorie  is  referenced  often  as  a  gram-calorie. 
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cross-section  area  A  being  illuminated  at  each  instant  of  time.  This 
arises  not  only  because  of  the  earth's  revolution  about  the  sun,  but 
also  because  of  satellite  tumbling.  For  the  purposes  of  the  present 
report,  only  spherical  satellites  will  be  considered. 

B.  DIRECT  RADIATION  FORCES  ON  A  SPHERE 
1.  Specular  Reflection 

Consider  a  spherical  satellite  of  radius  Rg  which  reflects  in  a  pure 
specular  manner  all  light  it  does  not  absorb,  having  a  uniform  re¬ 
flectivity  coefficient  kg.  The  incident  radiation  pressure  jd,  having 

A 

magnitude  p  and  direction  i  ,  can  be  taken  to  approach  the  satellite 

A  A°A  A  A 

in  the  coordinate  triad  (i, j,k)  shown  in  Figure  5  such  that  igg  - i. 

The  angle  is  defined  as  the  angle  that  £  makes  with  any  differential 

A 

element  of  surface  area  dA,  which  has  magnitude  cLA  and  outward  normal  i^. 

By  the  definition  of  specular  reflection,  the  luminous  flux  leaving  the 

surface  element  also  makes  an  angle  with  the  normal  in  accordance 

n 


with  Snell's  Law,  and  has  an  irradiance  I  =  k  I  . 

r  s  s 


FIGURE  5 .  SPECULAR  REFLECTION  FROM  A 
SPHERICAL  SATELLITE 
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In  spherical  coordinates. 


dA  =  R  2  sin  ©  d©  d^>, 
s 

"i^  =  (sin  ©  cos  (p)  ^  +  (sin  ©  sin  )  j  +  (cos©)  ic, 

^  ^  =  (2  sin2  ©  cos2^  -1)  'i'  +  2  (sin2  ©  sin^  cos(/> ) 

+2  (sin©  cos©  cosi$)1c. 

T  (1.7) 

/N  /\ 

The  last  unit  vector  can  be  computed  by  noting  that  i,  i^,  and 
ire^  ar®  necessarily  co-planar  and  by  working  with  the  simple 
sketch  in  Figure  6,  below. 


FIGURE  6.  GEOMETRY  FOR  DIRECTION  COSINES  OF  i 

reix, 


Once  the  relations  in  (A. 7)  are  used  to  obtain 

ire^  •  dA  =  Rg2  sin2  ©  cosi^d©  dp, 
equation  (1.4)  can  be  evaluated  over  the  illuminated  surface, 
the  hemisphere  defined  by  0  £  ©  £  n  and  -(tt/2).£.  <p  ^  (tt/ 2) , 
to  yield 

F,  =  -d  (I  /c)  R  2 
-d  v  s'  '  a 


tt  tt/2 

n 


Jsin2©  cos^> -k  (2  sin ^©  cos-'p  - 


0  -tt/2 

-sin2  ©  cosp  ")]  dp  d©  +  $  Odd((^)  4c  Odd  (©), 
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where  Odd  (<j>  )  is  the  integral  of  an  odd  function  of  (j>  ,  sin  <f>  • 
cos  <p  ,  over  the  range  -tt/2  to  tt/2,  and  Odd  (0)  is  the  integral 
of  a  function  of  0,  sirr^O  cos©,  which  is  odd  over  the  range  0  to 
tt.  Hence 


Odd  (  <f>  )  -  Odd  (0)  -  0, 


and  the  net  force  is  in  the  -i  direction.  Moreover;  the  factor 
multiplying  k  in  the  first  integral  is  clearly  odd  over  the  inter¬ 
val  -tt/2  £  (p  ^  tt/2.  Hence  if  F^  is  divided  by  the  mass  m,  and  the 
final  integration  performed. 


M 

•^■direct 


i8.  <Vc>  (V-), 


(1.8) 


2 

where  A  =nl  .  Note  that  the  constant  K  introduced  in  (1.5) 

8 

is  1.0,  independent  of  the  reflectivity  k,  for  this  pure  specular 


case  of  a  spherical  satellite. 


2.  Diffuse  Reflection 

Consider  a  spherical  satellite  having  uniform  reflectivity  kd 
which  reflects  in  a  pure  diffuse  manner  all  light  it  does  not  ab¬ 
sorb.  The  net  differential  forces  dF  are  normal  to  each  incremental 

-refl 

illuminated  area  dA,  i.e.,  in  the  direction  d^,  and  the  reflection 
obeys  Lambert's  Law.  The  magnitude  of  dF^^  can  be  deter¬ 
mined  by  centering  a  hemispherical  "cake  cover"  of  unity  radius 
over  dA,  as  suggested  in  Figure  7,  below,  and  by  noting  that  the 
total  luminous  flux  (power)  P  passing  through  the  cake  cover  is 
all  that  dA  emits:  it  is  k^  times  the  flux  incident  on  dA.  That  is, 

if  I  ( y )  is  defined  as  the  irradiance  in  the  direction  of  dS  due  to 
reflection  from  dA, 
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* 


FIGURE  7.  DETERMINATION  OF  NET  PRESSURE  NORMAL  TO 
INCREMENTAL  SURFACE  AREA  DUE  TO  DIFFUSE  REFLECTION 

D  -  M,  d..  •  <*±>  -  f  Uf)ds. 

cake  cover  ' 


U^) 


Denoting  the  irradiance  normal  to  dA  by  1^,  Lambert's  Law  provides 
for  diffuse  reflection 


I  (  y  )  -  IA  cos  'f  .  (1.10) 

Putting  this  into  the  integral  in  (1.9),  noting  that  dS  =  (sin^  JdV'd 6  , 

and  integrating  over  the  range  0  J  2rr,  we  find  that 

IA  - -(l/n)  kdIs  (is3  •  dA).  (1.11) 

The  pressure  on  the  cake  cover  due  to  the  light  reflected 

from  dA  is  clearly 

Erefl  =  (l/o)  1  <f)%>  (1.12) 


of  which  only  the  component  normal  to  dA  produces  a  reaction  force 
on  dA: 


dF  _ 
-refl 


-i. 


(i* 


J 


cake  cover 


«). 
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(1.13) 


From  (1.10),  (I.ll),  and  (1.12),  this  becomes 

d£refl  (a/3)  (kdVc)  (^ss  ‘  d±)* 

Hence,  in  (1.3), 

ir  =  (2/3)  kdis  des  •  tk).  d.H) 

Using  (1.7)  and  the  convenience  1  =-k,  which  is  the  unit 

ss  7 

vector  in  the  minus  vertical  direction,  (1.13)  can  be  integrated 


over  0  £  0  £  n/2  and  0  ■£  0  £  2rr  to  provide 

2nr  tt/2 


-refl 


-  iss  (2/3)  (kdIs/c)R^  J  J  cos2  0  sin  6  d  0  d  0 


0  ~  0 


‘is  U/9)  (kdIs/c)  nR2 


(1.15) 


Added  to  =  F^,  which  is  still  valid,  and  then  dividing  by  m, 


so  that  in  this  pure  diffuse  case  the  K  introduced  in  (1.5)  is 


(1.16) 


dependent  upon  the  reflectivity  constant  and  assumes  a  maximum  of  1.44. 


C.  SHADOW  REGIONS 


The  above  equations  depend  upon  the  direct  solar  radiation  pressure, 
which  is  directly  proportional  to  the  free-space  solar  radiation  pressure 
in  full  sunlight,  some  fraction  thereof  in  the  penumbra,  and  zero  in 
the  umbra. 

In  this  section  it  will  be  shown  that  the  penumbral  region  can  be 
neglected  and  the  umbral  regions  boundary  can  be  expressed  as  the  simple 
relationship 

r  |sin  a  |  =  R,  (rr/2)  £  a  £  3tt/2 
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where  r  is  the  geocentric  distance  to  the  satellite,  a  is  the  geocentric 
angle  between  the  sun  and  the  satellite  and  R  is  the  radius  of  the  earth  . 

The  use  of  this  approximation  results  in  errors  of  less  than  0.3  percent 
for  orbits  within  2000  miles  of  the  earth. 

Specifically,  it  suffices  to  assume  that  the  sun  and  earth  are 
perfect  spheres  and  the  satellite  a  point.  The  geometry  then  appears 
as  in  Figure  7 ,  below, 

R  „  =  radius  of  the  sun,  432  x  103  miles , 
sun 

R  =  radius  of  the  earth,  3*964  x  10^  miles, 

r  =  distance  from  sun  to  earth,  having  a  mean  R  of  93  x  10^  miles, 

6  S  ^  6S 

a  perihelion  of  91.6  x  10  miles,  and  an  aphelion  of 
94*7  x  10^  miles, 

r  =  radial  height  of  the  umbral  cone  measured  from  earth, 
c 

Penumbral 


FIGURE  8 


rc '  =  radial  distance  to  the  apex  of  the  penumbral  cone  as 

measured  from  earth. 
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By  similar  triangles, 

r  =  R  r  /  ( R  -R ) 

c  esM  sun 

=  a  mean  of  860,000  miles  and  850,000  ~  878,000  miles. 

Let  the  satellite  be  restricted  to  within  2000  miles  of  the  earth. 

Consider  the  vertical  line  marked  h£  in  Figure  8  to  cross  the 

sun-earth  line  at  a  distance  of  2000  miles  from  the  earth's  surface. 

The  ratio  of  he  to  R  is  clearly 

h/R  =  (r  -R-2000)/r 
c  c 

-  0.993. 

Therefore,  for  all  orbits  within  2000  miles  of  the  earth,  the  umbral 
region  can  be  represented  by  boundaries  parallel  to  the  earth-sun  line 
with  an  error  of  less  than  0.7  percent. 

For  the  penumbral  region,  similar  calculations  yield 

r'c  ^  845 ,000  miles , 

and 

h'c  =  (r'c  +  R  +  2000)/r'c  6  1.007. 

Hence,  for  a  satellite  orbit  in  an  earth-sun  plane,  the  penumbra  will 
be  less  than  0.7  percent  larger  than  a  cylindrical  shadow. 

The  use  of  a  cylindrical  umbral  region  with  no  penumbral  region 
at  all  involves  the  shrinking  of  the  region  of  partial  shadow  and 
expanding  the  region  of  total  shadow.  Except  in  certain  somewhat 
pathological  cases,  these  changes  will  tend  to  cancel  each  other  for 
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near-earth  orbits.  A  rule-of-thumb  estimate  indicates  that  the  maximum 
error  in  solar  pressure  due  to  the  approximation  will  be  on  the  order 
of  0.3  percent  on  any  single  pass,  with  the  effect  averaging  to  trivial 
proportions  on  multiple  passes.  In  any  event,  the  errors  due  to  the 
assumption  of  only  a  cylindrical  umbral  region  will  be  negligible  in 
comparison  with  the  uncertainties  in  solar  intensity,  in  refraction  of 
the  shadow  boundaries,  and  with  similar  ignorances  of  the  physical  situation. 

To  compute  when  the  satellite  is  within  these  simplified 
boundaries  is  just  a  matter  of  analyzing  the  geometry  of  Figure  9: 
the  satellite  is  in  complete  shadow  if 

|  sin  a  |  i  R/r,  (tt/2)  £a  £3^/2,  (1*17) 

and  in  full  sunlight  otherwise.  An  alternative  expression  avoids 

A 


FIGURE  9 

the  computation  of  a  from  inverse  trigonometric  functions .  Here 
(AJ.7)  appears  as 

complete  shadow:  J  sin  aj  —  R/r  and  cos  a  i.  0, 

full  sunlight:  otherwise.  (I. 17') 
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The  cos  a  can  be  determined  in  terns  of  the  ecliptic  coordinates  of  the 
»un  and  the  instantaneous  orbital  elements: 

y/r 
z/r 

where  =  celestial  longitude  of  the  sun  (measured 

in  the  ecliptic  plane  from  the  vernal  equinox) , 
i#  “  obliquity. 

Sin  a  then  follows  from 

sin  a  =  (1  -  cos  a) 

a/ 

Musen*®  obtains  expressions  for  the  long-term  effects  of  solar 
radiation  pressure  by  neglecting  shorter-term  effects  such  as  those  due 
to  the  shadow.  By  examining  long-term  effects  only,  he  is  able  to  estimate 
such  things  as  satellite  lifetimes.  However,  these  approximate  methods 
are  not  applicable  to  observations  on  only  a  few  periods. 

D.  REFLECTION  AND  RE-RADIATION 

Terrestial  radiation  or  "earth  shine"  pressure  also  exists .  Of  the 
total  insolation,  an  average  of  36  percent  is  reflected  or  back-scattered 
and  64  percent  is  absorbed  and  re-radiated  thermally."^  Reflection  and 
back-scattering  varies  between  15  percent  for  a  clear  sky  and  55  percent 
for  an  overcast  sky.  That  which  is  absorbed  and  re-radiated  as  heat  is 
primarily  counter  to  the  central-force  acceleration  of  gravity  and  virtually 
uniform  over  the  surface  of  the  earth.  Since  this  cannot  exceed  10-^g 


=  (cos  Xg)  x/r 
+(sin  \e)(cos  ig) 
+(sin  \ )(sin  i  ) 

•  Qf 
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for  the  satellite  and  orbit  types  specified  under  this  contract,  where  g 
is  the  acceleration  of  gravity,  it  is  doubtful  that  this  radial  component 
could  ever  be  identifiable  in  a  satellite  orbit. 

Hence  thermal  re-radiation  will  be  neglected  and  only  reflection  and 
back-scattering  will  be  examined.  These  two  reflections  have  the  alter¬ 
native  nomenclature  "specular"  (mirror-like)  reflection  and  "diffuse” 
reflection.  Of  the  two,  diffuse  reflection  by  far  comprises  the  major 
form  of  "earth  shine." 

1 .  Earth  Shine  Theory 

Consider  a  spherical  earth  illuminated  by  the  sun.  The  fireball 
on  the  earth’s  surface  can  be  neglected  because  of  its  small  contribution 
to  the  total  reflected  irradiance  at  satellite  altitudes  (viz,  pictures 
of  the  sunlit  earth  published  by  NASA  which  were  taken  by  its  Applications 
Technology  Satellite  (ATS)).  Since  the  fireball  is  the  manifestation 
of  the  specular  part  of  the  reflection,  its  negligibility  permits  the 
assumption  that  the  earth  is  a  perfect  diffuse  reflector. 

Much  of  the  analysis  performed  earlier  in  this  appendix  therefore 
applies.  Instead  of  applying  Lambert’s  Law  to  the  satellite,  however, 
we  now  apply  it  to  the  earth. 

Examine  Figure  10,  below,  where  for  convenience  we  have  positioned 

A  /%  A 

the  satellite  along  the  k  unit  vector  in  the  coordinate  triad  (i,  j,  k), 

A  a  A 

and  the  sun  has  been  put  along  the  i  ,  unit  vector  in  the  (i,  k)  plane. 

s  s 

We  maintain  the  convention  that  the  earth-sun  angle  a  lies  between  0 
and  nr.  The  unit  vector  i^  is  the  outward  normal  from  any  infinitesimal 
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'  A  ’  d  - 


A 

•  j 


NIGHT  SIDE  OF  EARTH 


FIGURE  10.  GEOMETRY  FOR  DIFFUSE 
EARTH  SHINE 
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element  of  area  dA  on  the  earth's  surface,  L  is  distance  between  the 
element  dA  and  the  satellite,  x^  is  the  element-to-satellite  unit  vector , 

A  A 

and  if  is  the  interior  angle  between  i  and  i^  . 

Cl early 9 

y\  A 

i  =  k  , 

r  ’ 

A  .  A  /  A 

i  =  -(sin  a)  i  -  (cos  ajk , 
s  s 

^i.  =  (sin  0  cos  $)i  +  (sin  0  sin  0)j  +  (cos  0)k, 

dA  =  R2  sin  0  d0  d0.  (1.19) 

With  1^  defined  as  the  component  of  irradiance  normal  to  dA  at  dA, 
we  employ  an  equation  much  like  (I.ll): 

IA  =  -  (1/n)  qlg  (iss  *  dA),  (1.20) 

where  q  =  earth's  albedo,  or  (diffuse)  reflectivity  coefficient.  From 

A 

Lambert's  Law,  the  component  in  the  direction  will  be  1^  cos  ^ 

At  the  satellite,  a  distance  L  away,  the  irradiance  l(  )  due  to  dA 
will  therefore  be 

K  ^  )  =  (ia/l2)  cos  /  •  (1.21) 

The  distance  L  can  be  computed  from  the  simple  sketch  given  in 
Figure  11,  below,  where,  by  inspection, 

a  =  r  cos  0  -R, 
b  =  r  sin  0  , 

L  =  J  a2  +  b2  =  r 


Jl  +  X2  -2X  cos  0  , 

(1.22) 


where  X  =  R/r. 
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FIGURE  11.  SKETCH  FOR  COMPUTING 
SURFACE-TO-SATELLITE  DISTANCE  L. 


Now  we 
the  vector 


develop  by  subtracting  the  geocentric  vector  to  dA  from 
r  =  rk  and  dividing  by  the  distance  L: 


As 


_  -X  (sin  0  cos  0)*i  -  X  (sin  Q  sin  ^)  .1  +  (l  -  X  cos  Q)  *k 


I 


+  X  -2 X  cos  0 


(1.23) 


Then,  from  (1.19), 


CCS 


cos  0  -X _ 

J  1  +  X*  -  2X  cos  0 


(1.24) 


Because  of  the  symmetry  about  the  (i,  "k)  plane,  the  solar  pressure  at 


the  satellite's  location  due  to  dA  has  its  component  out  of  that  plane 
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cancelled  by  another  area  element  on  the  other  side  of  the  plane.  The 
in-plane  components  are,  of  course,  doubled  in  magnitude  by  that  other 
area  element.  Hence,  restricting  0  to  be  between  0  and  rr,  we  find  that 
the  incremental  pressure  at  the  satellite  due  to  the  reflected  light  from 
dA  is 


dE  (8,  a)  -  2  [l(  ^  )/c]  (i'As)1;k, 


0  0  ^  TT 


(1.25) 


A  A 


where  (i^g)^  ^  =  the  projection  of  i^g  on  the  (i,k)  plane. 

For  simplicity,  now,  let  us  normalize  the  solar  pressure  by  removing 

the  environmentally  dependent  quantity  (ql  /c).  We  define  the  residual 

s 


expression  as  the  differential  illumination  factor 
dC(0,  0,  a)  =  djg(0,  0,  a)/  (qlg/c). 


(1.26) 


From  (1.25),  we  find  that 

dC(0  0  a)  =  (2/n)^c°3  6  zhlLsi n  9)(cos  a  cos  0  +  sin  a  sin  0  cos  0) 

(1  +  x2  _  2  X  cos  0)2 

•  £-X(sin  0  cos  0)  i  +  (l  -  X  cos  0)  dO  d0, 


which  we  can  express  more  compactly  as 

dC(O,0,a)  =  (2 !X2/tt)  [cos  a  (I^i+I^)  +  sin  a  (I^i+I^k)JdO  d0, 

(1.27) 

where 

r\ 

t  _  X  cos  0  sin^  Q  cos  0  (cos  Q  -X) 

^  (X2  +  1-  2X  cos  0)2 

j  _  cos  0  sin  0  (l  -  X  cos  0)  (cos  0  -X) 

2  (X2  +  1  -  2X  cos  0)2 
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X3  =  - 


X  sir?  0  cos2  0  (cos  0  -X) 

(X2  +  1  -  2X  cos  e)2 


j  _  sin2  Q  cos  0  (1-*  cos  Q)  (cos  6  -X) 

4  (X2  +  1  -  2X  cos  0)2 

The  total  earth-shine  illumination  factor  is  the  integral  of  (1.27) 
taken  over  that  part  of  the  earth's  sunlit  surface  that  is  visible  to  the 
satellite: 


C(a)  = 


Q1 

f 


.0, 


dC(O,0,a). 


(1.28) 


6  =  e  0  =  o 

O  r 


=  cos  (R/r)  =  cos  ^X , 


The  upper  limit  on  0  is  always 

(1.29) 

as  is  quickly  discernible  from  Figure  11.  The  lower  limit  on  0  and  the 
upper  limit  on  0  depend  on  the  satellite-earth-sun  geometry:  i.e.,  on 
a  and  X.  The  situation  decomposes  into  four  distinct  cases. 

Case  A.  The  Satellite  Sees  an  Earth  which  Is  Entirely  Sunlit 
(O<a.<n/2-0^). 

This  is  the  only  case  that  can  be  integrated  analytically.  Figure  10 
indicates  that  0q  =  0,  0^  =  n.  Hence,  from  (1.28), 

C(a )  =  X  (k  cos  a  +  i  J .  sin  a)  ,  (I. 30) 


where 


0, 


=  (  2X/tt)  J  J  I2d0  d0  =  -(l/2)(B1+B2B3+kB4), 
0  =  0  0  =  0 
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rr 


/ 

0  =  0 


I3  d0  d©  =  (1/4)  \ 

+<3k2+2kX-l)B3 


-(X+2k)(l-X)  + 
-(k+X)B4]  , 


Bj_  =  (l-X2)/2 
B  =  1-k2 

B3  =  log  I,  X  +  k  J 

n  _  n  ML _ i_ 

4  2  L  1  +  k  X+k 

k  =  (X2+l)/(-2X)  . 


Case  B.  The  Satellite  Sees  a  Dark  Region  Which  Covers  Less  Than  Half 
the  Visible  Earth  (rr/ 2  -  ©^  <  a  <  rr/2)  . 

We  approximate  the  sun  as  a  point  source  of  light  at  infinity.  Then 
the  terminator  (the  line  on  the  surface  of  the  earth  that  separates  night 
from  day)  is  a  great  circle.  The  terminator  intersects  the  horizon 
circle  of  the  satellite  at  the  two  points  A  and  B  shown  in  Figure  12, 
below.  If  a  is  less  than  90°,  i.e.,  if  less  than  half  the  visible  cap 
is  dark,  as  shown,  then  there  will  be  a  circle  (circle  0  in  the  figure) 
which  will  be  concentric  with  the  horizon  circle  and  tangent  to  the 
terminator  at  point  C  within  which  0  can  range  freely  from  0  to  ir. 

The  circle  is  at  0  =  rr/2  -  a. 

For  0  in  the  range  ©^  >  0  >rr/2  -  a,  the  range  on  0  will  be  a  constrained 
function  of  a  and  0,  having  upper  limit 

0j.  =  cos-^(-cot  a  cot  0).  (l.3l) 
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A  A 


CIRCLE  0 


HORIZON  AS  SEEN 
FROM  SATELLITE 


TERMINATOR 


FIGURE  12. 
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This  derives  from  Figure  13,  below,  wherein  the  following  relations 
apply: 

c  =  R  sin  0  sin  0, 


=  R  / 


1  -  c 


T 


R  J  1  -  sin^  0  sin2  0,' 


1  -  sin^  0  sin^  0  cos  a. 
The  length  e  can  also  be  obtained  as 
e  =  -  R  sin  0  cos  0. 

Equating  the  two  versions  of  e,  we  obtain  (I. 31). 


Now  (1.28)  becomes 


tt/  2- a  tt 


+  cos 


I' 


•/  /  (i^  +  kl2)  d0dO  + 


tt/  2-a 


C(a)  =  (2k<i/TT)|  k  cos  a  f  f  I^d0dO  +  i  sin  a  |  y 

*4=0  >0  cio  $=0 

*i  it 


I^d0dO+ 


n/ 2-a  0 


it 


+  sin 


taB  /  / 

n/2-a  7  0 


(il3  +  tcl4)  d0dO 


] 


(1.32) 


Case  C.  The  Satellite  Sees  a  Dark  Region  Which  Covers  More  Than 
Half  the  Visible  Earth  (tt/2  <  a  <  tt/2  +  0  ). 


We  refer  again  to  Figure  12,  except  that  now  the  terminator  lies  on 

A  A 

the  other  side  of  the  k  axis  (overhanging  the  positive  i  axis).  No 
Circle  0  now  exists  within  which  0  has  an  unconstrained  range  between  0 
and  tt.  Hence 
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FIGURE  13. 


C(a)  «  (2X2/n) 


(i  IL  +  kl2)  d0dO 


+  sin  a 


.  /* 

(il. 


,  kl4) 


d0dO 


] 


(1.33) 


Case  D,  The  Satellite  Sees  a  Completely 
Dark  Earth  (n/2 +©4  a  <  n). 
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No  earth  shine,  hence  C(a)  *  0.  Note,  also,  that  this  corresponds 
to  the  satellite's  being  in  the  umbral  region  specified  by  the  first 
line  of  (I.iy ).  Hence  there  is  no  solar  pressure  whatsoever. 

2.  Earth  Shine  Calculations 

The  evaluation  of  (I. 31)  and  the  numerical  integration  of  (1.32) 

and  (1.33)  yield  the  magnitude  and  angle  curves  of  Figure  14,  below. 

The  vertical  slash  on  the  right  end  of  each  curve  denotes  the  value 

a  =  rr/2  +  0^,  where  the  satellite  begins  to  see  a  wholly  dark  earth. 

Note  that  in  the  vicinity  of  a  equals  80°  or  90°,  the  magnitude  curves 

cross.  That  is,  the  illumination  for  large-enough  sun-satellite  angles 

is  not  a  monotone  decreasing  function  of  altitude.  This  contradicts 
29 

Dennison's  interpretation  of  his  previously  published  results,  which 
were  somewhat  more  complicated  than  those  of  Figure  14  and  provided  no 
angle  information. 

That  our  results  are  qualitatively  correct  is  clear  if  we  consider 
the  following  thought  experiment.  Put  a  satellite  just  off  the  surface 
of  the  earth  at,  say,  a  =  91°.  It  is  in  complete  dark,  and  hence 
C(a)  =  0.  Somewhat  higher,  at  a  100  km  altitude,  the  satellite  sees  a 
good  patch  of  sunlit  earth,  and  so  C( <1)7^0.  But  most  of  what  it  can  see 
is  dim  due  to  the  glancing  incidence  of  the  sun's  rays.  As  it  goes  yet 
higher,  it  sees  more  and  more  of  the  directly  lit  earth,  and  so  | [C(a) | | 
keeps  increasing  for  a  while.  Eventually,  however,  the  (R/r)^  dependence 
of  the  received  illumination  begins  to  take  effect,  and  | |C(a) | |  begins 
to  drop  with  further  increases  in  altitude.  Hence  the  illumination  for 
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FIGURE  14.  MAGNITUDE  AND  ANGLE  PLOTS  FOR  EARTH-SHINE 
ILLUMINATION  FACTOR  AT  SATELLITE. 
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large-enough  a  cannot  be  monotonic  with  altitude. 

For  purposes  of  digital  computation,  it  is  convenient  to  represent  the 
curves  of  Figure  14  in  an  approximate  analytic  form.  Since  the  direct 

A 

solar  radiation  impinges  in  the  direction  i  ,  the  approximation  will  be 

s  s 

A  A 

resolved  into  i  and  i  components ,  as  follows . 

S  S 


f  X  _  Z'  ^ 

C(aj  ~  q,  i  +  q_  i 
—  XL  r  ^2  ss- 


(1.34) 


where 


<h  = 


r X  (a^  +  a^)  cos  a, 
X  a 

r 

0 


M  1^/2-  ©1; 

,  tt/2-O-j^  ±  |  a  |  ^  n/2  +  0^; 

,  elsewhere; 


X  a. 


2  >  I  CL  I  <  tt/2  -  0X  ; 

x  ass  *  TT/2-0l  <  lal  £  n/2  +  ei  5 

(  0  ,  elsewhere  ; 


a1  -  (1/3) (-.0417  +  .5431  \), 

a ^  -  (1/3 )(. 0444-3. 17(X-. 77 )3+.0045(X-. 77 )sin  [l4.3(^-.77)rr]  ), 

ss  =  (3^2)  [l  +  s  -  s  esyr  -  e-T  y(  2  +  sy)]  , 

ar  ~  f  ass  +  (ai/2)  ls+1  -  s(l+sy)d]]  cos  a  + 

+  (1/6)  a)3  +  (X  -  sin  a)3]  _  il^£l!i;Jsln  a> 

l  (1  +  X2  -  2X  sin  a)3/2  X  J 

T  -  -4  +  9.3  ^  , 

y  =  (a  -  n/2)/01  , 


r- 1  ,  7  >  0  ; 

s  =l  1  ,  y<  0  ; 

d  =  3.7  +  59  (X-.77)2. 
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Figures  15  and  16,  below,  provide  a  comparison  of  the  approximate  and 

A 

exact  values  as  decomposed  into  the  k  and  -  i  (radial  and  cross-radial) 
components 


C^(a)  =  cos  a, 

-Ci(a)  =  q,  sin  a. 

E.  FINAL  RESULTS 

A  combination  of  all  the  details  for  a  spherical  satellite  yields 
finally 


r  ,  =  r,.  .  +  (ql  /c)C(a), 

-solar  -direct  s  —  * 

-  U3  I  (x,  t),  (1.35) 

where 

=  (1  +  4  kd/a)( Inom/c)(A/m)  is  a  solar 


"ballistic"  coefficient, 


=  qqx, 

P2  =  qq2’ 

all  other  terms  having  been  defined  in  the  body  of  this  appendix. 
Equation  (3.6),  in  the  main  text  of  the  report,  summarizes  all  these 
details  in  a  concise  way. 
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FIGURE  15.  APPROXIMATE  AND  EXACT  VALUES  OF  RADIAL 
EARTH-SHINE  ILLUMINATION  FACTOR 
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FIGURE  16.  APPROXIMATE  AND  EXACT  VALUES  OF 

CROSS-RADIAL  EARTH-SHINE  ILLUMINATION  FACTOR 
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F.  CLOSING  DISCUSSION 

Mainly  as  a  result  of  its  effect  on  the  orbit  of  Echo  I,  solar  radia¬ 
tion  is  regarded  as  the  most  significant  perturbative  effect  for  orbits 
whose  perigees  exceed  1000  km.  Large  variations  in  the  eccentricity  and 
geocentric  perigee  distance  for  such  orbits  are  almost  wholly  attributable 
to  the  effects  of  sunlight  pressure. 

Radiation  pressure  has  no  effect  on  the  period  if  the  orbit  is  circular. 

However,  if  the  orbit  is  non-circular  and  is  partly  in  shadow,  the  satellite 

can  enter  and  leave  the  shadow  region  at  different  distances  from  the 

sun,  resulting  in  a  net  gain  or  loss  of  energy  from  the  radiation  field. 

Even  if  the  orbit  does  not  pass  through  the  earth’s  shadow,  the  radiation 

pressure  has  the  effect  of  pushing  the  orbit  "sideways,"  so  that  its 

effect  on  the  perigee  does  not  vanish  even  for  circular  orbits. 

The  force  exerted  on  the  satellite  by  direct  solar  radiation  is  known 

to  within  about  2?  if  K(A/m)  is  known  perfectly.  Therefore,  the  effect 

of  the  radiation  is  highly  dependent  upon  the  accuracy  to  which  K,  A, 

and  m  are  known.  The  effect  of  thermal  terrestial  re-radiation  has  been 

neglected,  since  it  is  in  the  same  direction  as  that  of  the  principal 

gravitational  term  and  its  magnitude  is  negligible  in  comparison.  A 

16  18 

more  exact  analysis  is  available  in  Harvey  and  Fitz,  et.  al  .,  although 
this  accuracy  is  unnecessary. 
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APPENDIX  II 


ELECTRO: IAGNETIC  AND  ATMOSPHERIC  DRAG 

ON  SPHERICAL  SATELLITES 


DEFINITION  OF  SYMBOLS 


q  =  satellite  charge 

-19 

e  =  electron  charge  -1.6065X10  coulombs 

_n 

AQ  =  permeability  of  free  space  AtrXlO  (MKS  units) 

Jo 

=  permittivity  of  free  space  8.85X10  '  (MKS  units) 

M  =  earth's  magnetic  dipole  moment 

R  satellite  radius 
s 

r  =  satellite  distance  from  earth's  center 

Qg  =  satellite  potential 

v  =  satellite  velocity 
s 

B  =  magnetic  field  strength 

n  =  number  of  neutral  particles  per  unit  volume 


n^  =  number  of  ions  per  unit  volume 

m^  =  ion  mass 

mn  =  neutral  particle  mass 

b.  =  effective  satellite  radius 
i 
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A .  GENERAL 


There  are  basically  four  types  of  forces,  other  than  solar  radiation 
and  geopotential  forces ,  which  produce  effects  upon  earth  satellites . 

These  four  are  all  true  drag  forces.  In  order  to  determine  their  signifi¬ 
cance  it  is  necessary  to  have  a  reasonable  estimate  of  their  magnitude. 


A  satellite  accumulates  a  charge  q  while  moving  through  the  iono¬ 
sphere.  Since  the  satellite  moves  through  the  earth's  magnetic  field, 
it  experiences  a  force  given  by  the  basic  equation  governing  the  force 

19 

on  a  charged  particle  moving  in  a  magnetic  field.  According  to  Bechner  , 
the  magnitude  of  this  force  is 

q/i  Mv 

Fjjjgg  -  S  cos  (0-11.4°)  (II.1) 

2nr3 

for  a  polar  orbit  with  the  ascending  node  at  70.1°.  Here^Q  is  the 

permeability  of  free  space,  M  is  the  magnetic  dipole  moment  of  the 

earth,  v  is  satellite  velocity,  r  is  satellite  distance  from  the 
s 

earth's  center,  and  0  is  the  angle  between  the  radial  vector  to  the 
satellite  and  the  earth's  rotational  axis. 

In  order  to  perform  a  sample  computation  cos  (0  -  11.4°)  will  be 
assigned  its  maximum  value  of  1.  The  charge  q  is 

q  =  4tt  €0rsQs 
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where  £q  is  permittivity  of  free  space,  Rg  is  satellite  radius,  and 

Q  is  the  satellite  potential.  According  to  Brundin  in  Ref.  20,  Q 
s  s 

has  a  maximum  negative  potential  in  the  neighborhood  of  0.75V.  For 

computational  purposes  Qg  will  be  assigned  a  value  of  -1.0V. 

Assuming  a  spherical  satellite  of  radius  R  =  5.0m  at  an  altitude 

s 

of  200  km, 

q  =  4n( 8.85X10”^^) (5 )(l)  coulombs 

v  =  7.787X10^  m/sec 
s 


M  -K 

7 °  =  8.1X10  weber-m 

4n 


r  =  6.578X10°  m. 

The  pressure  is  Vmag/n^s^  equals  3  .li+X10~^n/m^. 


The  motion  of  a  conducting  satellite  in  the  earth's  magnetic  field 
will  cause  a  current  flow  and,  hence,  a  resultant  force.  This  force 
has  been  treated  by  Brundin^  and  found  to  have  a  magnitude  of 

~  eBR 

Find  ”  -eni  VRs  B(1"  - —  >  (II  -2) 

m.  v 
1  s 

below  the  hydrogen  region,  where  the  photoelectric  emission  has  little 
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effect  upon  the  effective  current.  At  higher  altitudes,  the  photoelectric 
emission  current  may  produce  forces  one  or  two  orders  of  magnitude  higher 
than  the  force  yielded  by  this  equation.  Because  of  the  difficulties 
in  predicting  photoemission  currents,  only  forces  at  lower  altitudes 
(  ft  800km)  will  be  considered.  For  a  more  detailed  discussion  of  photo¬ 
emission  current  see  Ref.  20. 

Here, 

e  =  electron  charge 

15 

n/K=  number  of  ions  per  cubic  meter  =  3*857X10  at  200km 

R  =5  meters 

3 

u  M 

B  =  magnetic  field  strength  =  /*o _ 

4nr3 

m^  “  ion  mass  =  2.5X10  kg  at  200km 

The  pressure  at  200km  is 


(1.6065)(3*857)(5)(8.1)X1011  \ A  $1.6065 )(8.1)(5)X10~4 
tt(6.578)3X1018  )\  (2.5)(7*787)(6.578)3XlO 


=  2.44X10  Un/m2 


D .  COULOMB  DRAG 

The  term  coulomb  drag  is  given  to  that  force  caused  by  incident 
ions  which  hit  the  satellite  because  of  the  satellite's  accumulated 
charge.  This  force  is  given  by 


F  .  =  nm.n.v^  (b<  -  R^) 
coul  1  i  s  i  s 


(11*3) 


2eQ  \l  . 

•  1  T 


where  b±  =  ^1+  2eW  s  j 


is  the  effective  radius  according  to  Ref.  20. 


Tn  ,,2 

m.  v 
l  s 


*  All  data  concerning  atmospheric  structure  was  obtained  from  Ref.  21, 
MODEL  10,  HOUR  0. 
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This  force  will  not  be  considered  independently,  as  it  is  a  component 
of  atmospheric  drag  and  will  be  treated  as  such. 

E.  ATMOSPHERIC  DRAG 

Atmospheric  drag  is  produced  by  the  collision  of  satellite  and  air 
particles.  It  therefore  includes  the  aforementioned  coulomb  drag.  It 
may  be  instructive  to  consider  the  atmospheric  drag  on  nonconducting 
as  well  as  conducting  satellites,  thereby  obtaining  some  appreciation 
of  the  coulomb  drag  effect. 


(a)  Atmospheric  drag  on  a  conducting  satellite 


where  is  drag  due  to  neutral  particles  and  is  due  to  incident  ions . 


If  nis  the  nuiiber  of  neutral  particles  encountered  per  unit  volume. 


F  =  nR^m  nv4 


(II. 4) 


n  s  n  s 


(II. 5) 


1  s 


Hence 
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and 


atm 


i”$C3)h 


(II.7) 


This  yields  a  reasonable  upper  bound  of 

Hi 


Fatm  -  (X  +  °-625  r  )  Fn 


(II. 8) 


since  m./m 
1  n 


2eQ 


0.5  and 


has  a  least  upper  bound  of  0.25  in  the 


_  ..2 
m.  v 
1  s 


region  being  considered. 

P  -  F  /rrR2 


n 


rr  s 

m  nv2 
n  s 


O 

where  v  -  GM  /r.  G  is  the  gravitational  constant  and  M  is  the  earth's 
sc  c 


mass.  (GM  =  3  .98866X101/*  m3/sec2) 


P  =  1.624  xlO  2  n/m2  at  200km 
n  ' 


and 


P  =  F/rrRj 
=  2. 26X10” 2  n/m2 

(b)  Atmospheric  drag  on  nonconducting  satellite 

The  force  equation  for  atmospheric  drag  on  a  nonconducting  satellite 
is  identical  to  that  for  a  conducting  satellite  (II. 7).  In  the  case  of 


the  nonconductor,  Q  is  zero  and,  therefore,  the  equation  reduces  to 

s 


W(l  +  °-5l)F„ 


(II. 9 
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Thus  the  ratio  of  nonconducting  drag  to  conducting  drag  is 


1  +  0.5 

n 

1  +  0.625  U± 
n 


(11.10) 


F.  RELATIVE  MAGNITUDES 

The  accompanying  table  and  figure  allow  comparison  of  the  four 
drag  effects .  In  addition,  the  solar  pressure  due  to  the  direct  radia¬ 
tion  of  the  sun  is  imposed  on  the  figure  for  the  cases  of  pure  specular 
reflection  and  pure  diffuse  total  reflection. 
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TABLE  1 


DRAG  PRESSURES  (n/m2) 
SPHERICAL  SATELLITES 


ALTITUDE  (km) 

AIRi 

AIR2 

INDUCED 

ELECTROMOTIVE 

CHARGE 

200 

2.26X10  "2 

2.1AX10"2 

2.A4X10-4 

3.14X10-12 

300 

2. 83X1 0"3 

2.5AX10"3 

5.29X10-5 

2.98X10-12 

500 

1.64X10"4 

1.36X10"^ 

4.81X10"6 

2. 68X1 0-12 

800 

6.91X10"6 

5.54X10-6 

2.26X10-7 

2.3XX10-12 

Air^  -  Atmospheric  drag  on  conducting  satellite 
Air^  -  Atmospheric  drag  on  nonconducting  satellite 

*  Helium  considered  to  be  charged  particles 
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DRAG  PRESSURE  (n/m‘ 


-i 

10 


102  km 

FIGURE  17.  Drag  Pressure  Versus  Altitude 
for  Spherical  Satellites 
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APPENDIX  III 


THE  EARTH »S  GEOPOTENTIAL  FIELD 


A.  INTRODUCTION 

The  potential  function  0(x)  for  the  geopotential  field  of  the  earth 
can  be  written  as  an  infinite  series  of  associated  Legendre  polynomials. 
In  truncated  form,  this  may  be  expressed  as 


where 


0M  =m/t  +  0o(x), 


N  n 

0  (x)  =  E  E  (C  U“  +  S  V®  ) 
0  n=2  m=0  m  n  ^  n 


(IH.l) 


(III. 2) 


is  the  potential  due  to  the  oblateness  (more  generally,  thp  asphericity) 
of  the  earth.  The  oblateness  acceleration  (x,J)  is  the  gradient  of 
with  respect  to  r: 


(III. 3; 


£°<i>  '  T  Jo  (C™»Br‘d  “n"  +  Sn-gr*d  O’ 

where  7  =  the  rotation  transformation  that  takes  the  earth- fixed 

geocentric  coordinate  system  (X1,  Y1,  Z1),  defined  as  the 
right-hand  system  with  X^"  at  Greenwich  and  along  the 
north-directed  polar  axis  at  time  t,  into  the  inertial 
geocentric  coordinate  system  (X,  Y,  Z)  defined  in  Figure  1, 
C  ,  ^nm  =  8eoP°tential  coefficients  from  (n,  m)  equal  to 

(2,  2)  upto  (N,  N)  whose  published  values  have  constant 
but  luiknown  biases  on  them, 

grad  =  column  (  ,  -^=-  , 

dx  dy  dzx 
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< 


u 


m 


n 


(n/r)  (R/r)n  Pnm 


(sin3 ) 


cos  mX 
sin  mX 


where  p  =  earth 1 s  gravitational  constant , 

r  =  the  geocentric  distance  to  the  satellite  at  time  t, 

R  =  earth’s  mean  equatorial  radius, 

Pnm(sin3)=  associated  Legendre  polynomials, 

3,  X  =  satellite  latitude  and  longitude,  respectively,  at  time  t. 


The  rotation  matrix  T  consists  of  a  precession  and  a  nutation 

necessary  to  align  the  polar  axis  at  time  t  with  the  polar  axis  at 

time  0  January  1,  and  then  a  rotation  around  the  polar  axis  to 

align  the  Greenwich  meridians.  (See  Reference  4.) 

27 

The  gradients  in  (ill. 3)  are  given  by^ 


grad  Un”  =  (l/R) 


lA’u'i-Ju'!1 

^  n  n+1  *-  n+1 


i,m  „  m-1 


Ur  mtl 


-  A  V  “  -  £v  “ 

^  n  n+1  n+1 


-  ( n-m+1 )  U 


m 

n+1 


5 


*We  use  the  definition 


^  W  ■  2nn! 


(l-x2)m/2  d' 


nt+n 


(x5  -l)n. 


dx 


nr+n 
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(III. 4) 


grad  Vnm  =  (l/R) 


i  Am  V  +:1 
n  n+1 


1  .m  T1  m-1  , 

2  A  U  , .  + 

n  n+1 


-(n-m+l)  V 


m 

n+1 


i  v. 


£  U 


m+1 

n+1 

m+1 

n+1 


where  A®  =  (n-m+l)  (n-m+2). 

Note  that  with  the  vector  J  defined  to  be  the  column  array  of 
geopotential  coefficients 

J  =  column  (C^,  C30,...  c22;  °31>  *  *  *  *  °33;  *  *  *  S22’  S31’  *  *  *  *  S33 ;  *  *  * J’ 

(III. 5) 


aF 

— O 

the  sensitivity  matrix  dJ  can  be  written  as 


dF 

=  T  [gradU2°,  gradU^0,...  gradU^2;  grad  U^1... 

gradll^;  ...  gradV,.2;  gradV^1,...,  gradV^;...] 


(III. 6) 


The  gradient  -  gradient  of  can  also  be  found 


27 


>%. 


dx^  dx^ 


aF 

-o 

ar 


i,j=l,2,3 


N  N 

=  T  E  E  (C  grad-grad  U™  ^ 

n=2  m=0  r,m  n 


+  grad-grad  vj)  T1 , 


grad-grad  = 


2R 


a“  (grad  U^)7  -  (grad  U^)T 
-a”  (grad  -  (grad  )T 

-2  (n-m+l)  (grad  l/^fl)T 


(HI. 7) 
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grad-grad  v”1  = 
n 


2R 


’An  Ci)T  -  («rad  ^i)T 

An  <«>*»  Cl)T  +  (S™1  C>T 
_  -2  (n-ntf-1)  (grad  V^fl)T 


Note  that  both  the  gradient  and  gradient-gradient  expressions 
require  the  negative-order  expression  for  efficient  computer  calcul¬ 
ation. 


((J-m, 

< 

=  (-l)m  _(n-m)_L  < 

U111  ^ 

\ 

n 

km> 

(nfm) ! 

-v111 

'  n  / 

(III. 8) 


Note,  for  purposes  of  numerical  checks,  that 


dF 

-o 

dr 


is  clearly 


symmetric  and,  moreover,  its  diagonal  elements  are  the  Laplacian  of 
0Q  and  hence  must  vanish  identically  in  free  space: 


2  2 

„2d  dr  ,  d*F 
\7  Y>  o  =  ox  +  oy 


2 

d*F 


oz 


a  0. 


dx 


dz 


( III. 9) 


3.  GEOPOTENTIAL  MODEL  COEFFICIENTS 


The  geopotential  model  to  be  used  will  be  the  essentially  eighth- 

23 

order  model  published  at  the  1966  COSPAR  by  Gaposchkin  and  reproduced 
by  V/ackernagel.  ^  It  consists  of  Kozai's  thirteen  zonal  coefficients 
C^q  =  q  =  -J.^,  plus  the  thirty-four  pairs  of  tesseral  co¬ 

efficients  (2,2),  (3,1)  to  (3,3),  •••,  (8,1)  to  (8,8),  plus  the  sixteen 
pairs  of  resonance  terms  (9,1),  (9,2),  (9,9),  (10,1)  to  (10,4),  (11,1), 
(12,1),  (13,12),  (13,13),  (14,1),  and  (15,12)  to  (15,14).  The  un¬ 
moralized  zonals  are  presented  in  Table  2  and  the  normalized  tesser- 
als  are  in  Table  3 • 
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The  term  "normalized"  applies  when  the  conventional  spherical 
harmonics  F^(x)  are  replaced  with  the  fully  normalized  spherical  harmonics 


Pjj(x)  =  !  ( 2n+l)  k/(n-Hn)l  F^(x)  =  N(n,m)F^(x) 


(III.10) 


where 


k  = 


1,  m  =  0; 

2,  m  ^  0. 


The  coefficients  C  ,  S  are  then  replaced  with  their  normalized  versions: 

run  nm 


“  cra/N(n-) 

-  Snr/N<n'm) 


(III.ll) 


wherever  the  Imx)  appear. 

The  geopotential-coefficient  covariance  matrix  Pj f  is  113  x  113, 
although  the  majority  of  its  elements  are  zero.  Tables  4  and  6  present 
the  standard  deviations  on  the  coefficients,  and  Tables  5  and  7  their 
correlation  matrix.  We  were  unable  to  fine  the  standard  deviations  for 

n 

the  (9,9)  terms  and  therefore  used  the  rule-of-thumb  value  0.2  x  10-' 

commonly  applied  to  these  normalized  coefficients.  This  pair  was  determined 

from  the  special  resonance  properties  of  MIDAS,  and  hence  are  uncorrelated 

23 

to  any  other  coefficients. 

The  coefficients  form  essentially  three  uncorrelated  groups:  the 
even  zonals,  the  odd  zonals,  and  the  tesserals.  The  correlation  matrix, 
therefore,  is  presented  in  three  parts,  one  for  each  of  these  groups. 
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Because  of  its  size,  the  tesseral  correlation  matrix  is  presented 
in  thirty-one  pages.  The  placement  of  these  pages  is  keyed  to  the 
correlation-matrix  map  in  Figure  18. 
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10 

J, 


2. 


14 


1082. 645  X  10  , 

-1. 649  X  10'°  , 

0. 646  X  10'6  , 

-0. 270  x  10"6  , 

-0. 054  X10'6  , 

-0. 357  X  10"6  , 

0.  179  X  10_b  . 

Table  2.  Unnormalized 


J3  =  -2. 546  X  10 

Jc  =  -0. 210  x  10 
b 

J7  =  -0. 333x10 
=-0.053  x  10 

J,  .  =  0.  302  X  10 

J13  =  -0. 114x10 

Zonal  Coefficients 
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a 

2 

1 

2 

3 

1 

2 

3 

4 

1 

2 

3 

4 

5 

1 

2 

3 

4 

5 

6 

1 

2 

3 

4 

5 

6 


n  m 


S  X  10" 


CX101 


S  X  10" 


C  X  106 

2.379 

1.936 

0.734 

0.561 

-0.572 

0.330 

0.851 

-0.053 

-0.079 

0.631 

-0.520 

-0.265 

0.156 

-0.047 

0.069 

-0.054 

-0.044 

-0.313 

-0.040 

0.197 

0.364 

0.250 

-0.152 

0.076 

-0.209 


-1.351 

0.266 

-0.538 

1.620 

-0.469 

0.661 

-0.190 

0.230 

-0.103 

-0.232 

0.007 

0.064 

-0.592 

-0.027 

-0.366 

0.031 

0.518 

0.458 

-0.155 

0.156 

0.163 

0.018 

-0.102 

0.054 

0.063 


7  7 

8  1 


8 

2 

8 

3 

8 

4 

8 

5 

8 

6 

8 

7 

8 

8 

9 

1 

9 

2 

9 

9 

10 

1 

10 

2 

10 

3 

10 

4 

11 

1 

12 

1 

12 

2 

13 

12 

13 

13 

14 

1 

15 

12 

15 

13 

15 

14 

0.055 

-0.075 

0.026 

-0.037 

-0.212 

-0.053 

-0.017 

-0.0087 

-0.248 

0.117 

-0.0040 

-0.065 

0.105 

-0.105 

-0.065 

-0.074 

-0.053 

-0.163 

-0.103 

-0.058 

-0.075 

-0.015 

-0.062 

-0.063 

0.0083 


0.096 

0.065 

0.039 

0.004 

-0.012 

0.118 

0.318 

0.031 

0.102 

0.012 

0.035 

0.0909 

-0.126 

-0.042 

0.030 

-0.111 

0.015 

-0.071 

-0.0051 

0.048 

0.010 

0.0053 

0.058 

-0.066 

-0.0201 


Table  3*  Normalized  Tesseral  Coefficients 
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Table  4.  Geopotential  Coefficient  Standard 
Deviations  (Unormalized)  -  Zonal  Coefficients 
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J2  J4 

J6 

J8 

J10 

J12 

J14 

J2 

1  *  00  —  .60 

,80 

-.89 

.79 

-.71 

.83 

J4 

1.00 

" « 86 

•  60 

-.85 

.91 

-.47 

J6 

1,00 

-.79 

.96 

CO 

CD 

• 

1 

.60 

J8 

1.00 

o 

CD 

• 

I 

.84 

-.84 

J10 

1.00 

— .  80 

.70 

J12 

1.00 

-.50 

J14 

1.00 

J3  J5 

J7 

J9 

Jll 

J13 

J3 

~1.00  -.93 

.98 

-.94 

.48 

— .  86 

J5 

1.00 

-.96 

.86 

-.69 

.75 

J7 

1.00 

-.92 

.57 

-.82 

J9 

1,00 

-.27 

.97 

Jll 

1.00 

-.12 

J13 

1.00 

Table  5 .  Geopotential  Coefficient  Correlation  Matrix  -  Zonal  Coefficients 
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FIGURE  18.  Geopotential  Coefficient  Correlation  Matrix  -  Tessera!  Coefficient  Reference  Chart  for  Table 
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TABLE  7.  GEOPOTENTIAL  COEFFICIENT  CORRELATION  MATRIX  -  TESSERAL  COEFFICIENTS  (PART  10  of  31) 
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TABLE  7.  GEOPOTENTIAL  COEFFICIENT  CORRELATION  MATRIX  -  TESSSRAL  COEFFICIENTS  (PART  11  of  31) 
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TABLE  7  .  GEOPOTENTIAL  COEFFICIENT  CORRELATION  MATRIX  -  TESSERAL  COEFFICIENTS  (PART  12  of  31) 
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TABLE  7.  GEOPOTENTIAL  COEFFICIENT  CORRELATION  MATRIX  -  TESSERAL  COEFFICIENTS  (PART  13  of  31) 
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TABLE  7.  GEOPOTENTIAL  COEFFICIENT  CORRELATION  MATRIX  -  TESSERAL  COEFFICIENTS  (PART  14  of  31 ) 
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TABLE  7-  GEOPOTTNTIAL  COEFFICIENT  CORRELATION  MATRIX  -  TFBSERAL  COEFFICIENTS  (PART  15  of  31) 
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TABLE  7.  GEOPOTENTIAL  COEFFICIENT  CORRELATION  MATRIX  -  TESSERAL  COEFFICIENTS  (PART  16  of  31) 
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TABLE  7.  GEO  POTENTIAL  COEFFICIENT  CORRELATION  MATRIX  -  TESSERAL  COEFFICIENTS  (PART  17  of  31) 
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TABLE  7.  GEOPOTENTIAL  COEFFICIENT  CORRELATION  MATRIX  -  TESSERAL  COEFFICIENTS  (PART  18  of  31) 
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TABLE  7.  GEO  POTENTIAL  COEFFICIENT  CORRELATION  MATRIX  -  TESSERAL  COEFFICIENTS  (PART  21  of  31) 
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TABLE  7.  GEOPOTENTIAL  COEFFICIENT  CORRELATION  MATRIX  -  TE5SERAL  COEFFICIENTS  (PART  22  of  31) 
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TABIZ  7 .  GEOPOTENTIAL  COEFFICIENT  CORRELATION  MATRIX  -  TESSERAL  COEFFICIENTS  (PART  26  of  31) 
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TABLE  7-  GEOPOTENTIAL  COEFFICIENT  CORRELATION  MATRIX  -  TESSERAL  COEFFICIENTS  (PART  29  of  31) 
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TABIE  7.  GEOPOT0ITIAL  COEFFICIENT  CORRELATION  MATRIX  -  TESSERAL  COEFFICIENTS  (PART  30  of  31) 


S12 » 2  S13, 12  S13» 13  514,1  S15*l?  Sl5,13  S15,14 


CM 

^Hl 

to 


in 


o 

H 

in 


OJ 

ro 

CNJ 

ro 

■H 

OJ 

fO 

OJ 

•H 

H 

•H 

OJ 

O' 

•> 

* 

- ^ 

* 

* 

o 

o 

o 

o 

OJ 

ro 

ro 

in 

in 

in 

O' 

O' 

•«« 

in 

in 

tO 

in 

in 

in 

l/) 

in 

to 

to 

to 

to 

to 

in 

in 

cm 

o 

OJ 

H 

ro 

o 

ro 

ro 

ro 

OJ 

O 

o 

G 

G 

O 

o 

G 

G 

o 

O 

o 

o 

O 

o 

O 

O 

1 

1 

i 

1 

1 

1 

c 

*-* 

cm 

ro 

OJ 

in 

o 

o 

CNJ 

CNJ 

G 

CNj 

G 

O 

o 

o 

o 

o 

o 

o 

o 

o 

o 

O' 

o 

o 

O 

1 

i 

i 

1 

f 

1 

— 1 

G 

CM 

•Ef 

ro 

O 

CM 

o 

ro 

H 

o 

G 

C 

G 

G 

C 

G 

G 

G 

C 

o 

o 

C 

o 

1 

1 

1 

r 

w— < 

1 

1 

H 

•— i 

G 

ro 

CM 

OJ 

o 

•— i 

CNJ 

«r-» 

•H 

c 

o 

O 

r- 

o 

O 

o 

o 

O 

o 

o 

o 

1 

1 

1 

1 

1 

1 

H 

o 

o 

ro 

CM 

O 

ro 

o 

O 

o 

G 

C 

G 

G 

G 

G 

o 

o 

o 

1 

1 

1 

1 

1 

1 

l 

CM 

o 

CM 

ro 

o 

CM 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

1 

1 

1 

1 

•H 

& 

O 

H 

O' 

ro 

o 

ro 

H 

o 

H 

G 

C 

o 

o 

o 

o 

O 

o 

1 

1 

1 

1 

1 

4P-4 

•H 

o 

O' 

ro 

o 

CC 

o 

o 

o 

r- 

i 

o 

1 

o 

o 

1 

o 

o 

4-1 

c 

in 

CM 

1^ 

O 

o 

G 

o 

O 

o 

G 

o 

• 

• 

• 

• 

• 

• 

• 

1 

i 

CM 

G 

ro 

sO 

CM 

C 

o 

O 

o 

O 

o 

# 

• 

• 

• 

• 

• 

#  H 


ro 

4k 

o 

in 


OJ  O  CM 
c  c  c 

•  •  • 


o  o 

G  O 

•  • 

I  H 


CM 

O' 

CM 

ro 

•H 

4H 

CM 

CM 

ro 

■H 

CM 

ro 

rH 

,14 

O 

o 

o 

o 

H 

CM 

CM 

IO 

ro 

* 

in 

in 

in 

O' 

O' 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

H 

in 

in 

in 

in 

in 

in 

in 

in 

in 

in 

in 

in 

in 

in 

in 

136 


TABLE  7-  GEOPOTLNTIAL  COEFFICIENT  CORRELATION  MATRIX  -  TESSERAL  COEFFICIENTS  (PART  31  of  31) 


APPENDIX  IV 


TRANSITION  MATRIX  $(t,t  )  IN  RECTANGULAR 
COORDINATES 


Repeating  equations  (2.2), 
r  =  r0f(AE)  +  ^  g(  A  E) , 

r  =<ir/dt  =  Fq  ft(  A  E)  +  Jq  gt(  A  E),  (IV.1) 

where  the  f  and  g  functions  are  defined  with  equations  (2.2).  From 
(2.6),  the  Keplerian  transition  matrix  is 


dr/  &£o  | 

br/  bio 

d£/  b  Iq  1 

bt/  dlo 

$(t,tQ)  = 


evaluated  on  the  two-body  orbit  (IV.1),  where  the  state  vector  is, 
of  course } 

r 
f 


(IV. 2) 


From  (IV.1),  the  four  3x3  submatrices  are 
dr/  b  Tq  =  fl  +  r^  d  f/  b  Tq)  +  Tq  (  bg/  d  Tq)  , 
br/  b  r^  =  i^(  b  f/  b  Tq)  +  gl  +  (  b  g/  b^)  , 
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df/  a Sg  =  ftl  +  Tq  (  dlt/  drj  +  i^(  dgtZ  d %)  , 

dt/  t>*0  -  £o(  dV  dk>)  +  St1  +  io  (  dh/  • 

where,  for  example,  the  row  vector 

(  df/  di^)  =  (  df/  dxQ1,  df/  dxQ2,  df/  dx^) 

Noting  that 


xi+3 


1  -  i  *  3, 


we  can  employ  for  the  sake  of  clarity  the  convention  that  the  x 
subscripts  will  never  exceed  3  and  derivatives  will  be  denoted  by 
(")  rather  than  by  the  subscripts  A  to  6.  Then 


1 

X2 

,  f  = 

*2 

x„ 

L  3j 

-*3- 

r  =  (x12  +  x22  +  x32)  , 


r0  =  (x01  +  ^2  +  X03) 


r  =  (ij2  +  x22  +  ±32)\ 


X01 

X01 

X02 

»  ^0 

*02 

_^3_ 

_x03_ 

.  _  , .  2  ,  .  2  ,  .  2\; 

r0  ~  (x01  +  ^2  +  X03 


(IV  .3) 
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In  these  terms,  we  can  write  the  semi-major  axis  as 


=  *vQ  /  (2 M  -Vo 


0*0 


the  inner  product 


d0  =  io  •  So  =  *01  *01  +  *02  *02  *  *03  *03 


will  also  be  useful. 


From  the  definitions  of  the  f  and  g  functions, 
ax0i -  r02  (  d  a/  dxoi} 


df/  = 


df/d*0i 


r3 

0 


d(AE) 

(1-cos  A E)-(a/rQ)sin  AE  dx^  ’ 


(  5  a/  diu. )  d(AE) 

- - -  (1  -  cos  A  E)-(a/rQ)  sin  AE  ‘ 


l0 


dX0i  ' 


(pa)*x  .  sin  AE  +  ar-x,.  (1-cos  A  E) 

b  8/  6  *bi  =  - - -  +  Vd  a/dx0i)+A5 


>^0 


5g/  d  iQi  =  (a/p)x0i  (1-cos  A  E)  +  A1  (  d  a/  & 

dxOi 


^ft^  b  *Oi 


_  _  d(r/a)  _  da  /\^ 

"  A-  [^i  +  2«0  3^-  +  a^7]  - 


cosAE 


a(AE) 

dx0i 


dVdiDi  ■  A2rot21 


d  (r/a) 


d  a 


)2 


bgt/ 

a  gt/  a  x0i 


2  ‘0L‘“  a  xrji  a  Xqi  rr 

a  (rA) 


a(AE) 


+  (r/a)  ]  -  ImI_  cos 

0  d*0i 


=  A,  — -v-A-AA— A  -(a/r)  sin  AE  ^(aI) 


3  d  *bi  . .  “  d*oi  ’ 

=  A  — “(a/r)  3in  AE  _ 

3  ^  <*01 


(IV  .4) 
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The  partials  that  appear  in  (IV. 4)  are  simply 


Sa/d^i 


*01  ^  +  s‘:o  > 

ro<:^  -  Vo2) 


ar'/dy;ji  = 


2^o  *01 

(^-Vo2)2 


d(r/a)  p  Xq^cos  AE  +  a'  rQxQj_san 


6x0i 


AE 


F  ar0 

+  A6  (6AE/dx0i), 


-  A^( da/ axQi ) 


d(r/a)  x(  sin  AE 
dx0i  (^a); 


where 

dAE 


dx0i 


-  A^(da/dxQi) 


i-  A^(dAE/dxQi)  , 


=  Ar 


*01  -  iffE1 

0  (pa)- 


*6i 


6AE 

Ox0i 


=  A 


A  (da/6io  ) .  ^ 

(/aa)- 


and  the  parameters  A,  through  A^  are 


Al  = 


p  r^sin  AE  +  2a'dQ(l-cos  AE) 


2pa 


,  _  (pa) ‘  sin  AE 

2  j*2*  3 

,x  r0 
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(a/r)^(l-cos  A  E), 


^^r^cos  AE  +  a^d^  sin 

2)iK2 


AE 


8 


(au)2  Tq  cos  AE  +  ad^  sin  AE 

~  , 

(l  -  ^/a)  sin  AE  +  d^cos  AE)/(ua)£  , 

(r^/a2)  sin  AE  +  d0^  ~  c°s  ^  _  (3/2)(u/a^)2  At, 

2  ( ua^ )  ^ 

_ 1 _ 

1  +  d^  (sin  AE)/(ua)£  -  (l-r^/a)  cos  AE 


•  o 

Clearly,  when  r^r^  approaches  2  )i,  the  elements  of  the  transition 
matrix  tend  to  be  ill-conditioned  because  of  a  and  its  partials.  This 
corresponds  to  a  near-parabolic  orbit.  Although  this  case  does  not  apply 
to  the  present  study,  it  was  considered  in  the  development  of  MINIVAR, 
and  hence  the  so-called  NASA  orbital-element  states  are  used  there.  Only 
one  of  the  states  reflects  a,  and  the  ill-conditionedness  of  the  corresponding 
transition  matrix  is  virtually  avoided  in  the  near- parabolic  case. 

To  go  from  the  transition  matrix  of  this  appendix  to  the  form  used 
in  MINIVAR  requires  only  a  point  transformation,  as  outlined  in  Reference  4. 

A  comparison  of  this  sort  was  programmed,  and  in  all  cases  the  two  transition 
matrices  agreed  to  within  round-off  tolerances.  Thus,  both  the  algebraic 
details  developed  here  and  the  MINIVAR  development  are  sustantiated. 
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APPENDIX  V 


COMPARISON  OF  MAXIMUM  LIKELIHOOD  AND 
MINIMUM-VARIANCE  ESTIMATION  OF  SPACE-VEHICLE  MASS 

Suppose  that  the  sequence  {x(^)j  k  =  0,1,...  of  real 
random  n-vectors  x(k)  is  governed  by  the  recursive  equation 

x( k+1 )  =  $(k+l,  k)  x(k)  f(k)  +  ag(k),  (V.l) 

where  (^)  (k+1,  k)  is  a  given  n  x  n  transition  matrix,  a  is  an 
unknown  system  parameter  (such  as  the  area-to-mass  ratio  of  the 
main  report)  which  will  be  regarded  as  a  real  random  variable  having 
a  priori  mean  a  and  variance  C~  ;  f(k)  and  g(k)  are  random  n-vectors 

with  means  ?(k)  and  g(k)  respectively,  and  x(0)  is  a  random  n- vector 
with  mean  x(0)  and  covariance  matrix  Pq.  In  addition,  suppose 
that  for  each  k  =  1,  2,  . ..,  N,  we  have  available  an  m  x  1 
observation  vector  z(k)  given  by 

z(k)  =  H(k)x(k)  +  v(k) ,  (V.2) 

where  H(k)  is  a  given  m  x  n  matrix  and  v(k)  is  a  random  m-vector 
of  observation  errors  having  mean  zero  and  covariance  matrix  R(k). 
Finally,  we  will  assume  that  all  of  the  vectors  x(0),  f(0),  f(l),..., 
g(0),  g(l ),...,  v(l),  v( 2) , . . . ,  are  pairwise  uncorrelated,  and  that 
the  random  parameter  a  is  strictly  independent  of  all  of  these  vectors. 
Using  the  available  observations,  we  wish  to  determine  an  optimal 
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(in  some  well-defined  sense)  a  posteriori  estimate  of  the  parameter  a. 

For  any  random  vector  y,  if  y  is  the  a  priori  (unconditional) 
mean  of  y,  we  will  write  2  =  Z~Z>  y  =  y  +  y  is  a  well-defined 

decomposition  of  y  into  a  deterministic  part  and  a  zero  mean 
random  part.  In  equation  (V.l),  we  then  define  e(k)  =  ?(k)  +  a|[(k), 
and  note  that  e(k)  is  orthogonal  to  x(0),  v(l),  v(2),...,  a,  and 
to  each  e(j)  for  j  k.  With  this  notation  we  have 

x(k+l)  =  $(k+l,k)x(k)  +  T(k)  +  afi(k)  +  e(k) .  (V.la) 

The  covariance  matrix  of  the  zero  mean  random  vector  e(k)  will  be 
denoted  by  Q(k). 

If  h^  and  h^  are  zero  mean  real  random  variables  with  finite 
variances,  we  define  the  scalar  product  (h^,  h  )  =  Eh^h^  (where 
E( . )  is  the  expectation  operator),  and  the  norm  It  h  II  =  (h,h)*-. 

Let  %  denote  the  Hilbert  space  which  is  the  closure  in  this  norm 
of  the  linear  manifold  generated  by  all  the  components  of  $?(0), 
e(l),  ®(2),...,  v(l) ,  v(2),...,  and  by  a.  Clearly,  for  each  k  =  1, 
2,...,  the  components  of  x(k)  and  z*(k)  are  elements  of  W  ,  since  each 
of  these  components  is  expressible  as  a  finite  linear  sum  of  elements 
of  the  generating  set.  For  each  positive  integer  N,  let  9??  (N)  be  the 
finite  subspace  of  V  which  is  spanned  by  the  components  of  #(l),..., 
f(N),  and  let  K~  be  the  orthogonal  projection  on  (N). 

If  h^  is  an  arbitrary  element  of  (N),  then 

lla-hjl  =  ||S  -  k/II  +  II  KjjU  -  II  ,  (V.3) 
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since  -  h^T  €  (N)  and  a  -  K^a  _L  (N)  by  the  orthogonal 

projection  theorem.  From  (E.3)  it  follows  that  is  the  unique 
element  of  (N)  which  minimizes  the  distance  from  a  to  ty(N)j  i.e., 
K^a  is  the  unique  (up  to  equivalence)  minimum- variance ,  linear 
estimate  of  a  with  respect  to  the  observations  J[(l),...,  £(N). 

If  the  components  of  zf(l) , . .  .z^N)  are  all  linearly  independent, 

then  the  projection  K^a  can  be  represented  as  follows.  Let  denote 

*T  p  a/T  /  \  /  \  t 

the  vector  whose  transpose  is  defined  as  =  |_z(l),...,z(N)J 

The  linear  independence  of  the  components  of  2^  means  that  the 

matrix  cov  (2^,)  =  e2^2^,  which  we  will  henceforth  call  U^,,  is  positive 

definite.  Setting  2^  =  ,  we  see  that  cov  (2^)  =  I.  Therefore, 

the  components  of  2j,  comprise  an  orthonormal  basis  for  Vi  (h), 

and  we  can  write 

-  0^)%  -  (es|h)t%%-  u.*) 

where  the  components  of  EfiC^are  the  Fourier  coefficients  of  a  with 
respect  to  the  orthonormalized  observations.  Since  &  =  a-a,  the 
minimum- variance,  linear  estimate  of  a  based  on  N  observation  points  is 
given  by 


a(N)  =  a  +  (EeC2JJ)TUN1(^J  -  \) . 


(V.5) 


-  144  - 


The  variance  of  the  estimation  error  may  be  computed  as  follows: 


E(a-&(N))2=  ii  s  -  i^anr 


=  <r 


-  (h»VTun1(e!1V  • 


(V.6) 


REMARK:  We  note  that  even  if  the  components  of  2^  are  not 

all  linearly  independent,  equations  (V.4),  (V.5)  and 
(V.6)  remain  valid  if  we  replace  U"1  by  U^,  the  generalized 
inverse  of  the  positive  semi-definite  matrix  Uj.  in  the 


10 


sense  of  Penrose  .  The  question  of  linear  independence, 
which  we  will  not  consider  in  this  appendix,  is  probably 
most  easily  discussed  using  the  notation  of  the  sequential 
estimation  procedure  of  Kalman,  described  below.  A 
sufficient,  but  by  no  means  necessary,  condition  for  the 
linear  independence  of  the  components  of  the  £(i)  s  is 
that  rank  (R(k))=  m  for  each  k. 

Let  us  suppose  for  a  moment  that  the  random  variable  £  and 
each  of  the  random  vectors  x(0),  e(0),  e(l),...,  v(l),  v(2),...,  has 
a  Gaussian  distribution.  In  this  case,  the  orthogonality  of  these 
random  quantities  implies  that  they  are  actually  strictly  independent 
of  one  another.  For  our  purposes,  the  important  fact  is  that  O' 
and  the  components  of  2^  will  then  have  a  joint  Gaussian  distribution. 
Hence,  it  is  quite  easy  to  obtain  the  maximum  likelihood  estimate  of 
a  given  2^.  For  this  and  later  computations,  we  will  need  the 
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following  result . 


LEMMA:  Let  the  symmetric,  positive  definite  matrix  M  be 
partitioned  as 

'a  b" 


where  A  is  pxp,  Bis  pxq,  and  C  is  q  x  q.  Then  the  matrices 
-IT  T  -1 

A,  C,  A  -  BC  B  and  C-B  A  B  are  each  positive  definite.  Further¬ 
more,  if  x  and  u  are  p-vectors  and  y  and  v  are  q- vectors,  then  the 
bilinear  form 


Q(x>  y,  U,  v)  =  [xT ,  yx] 

'a  b" 

-i 

u 

T 

— T 

O 

V 

admits  the  expansions 

(i)  Q(x,  y,  u,  v)  =  x^A'^u  +  (y  -  B^A'^x)^ C-B^A~^B)~^( v-B^A'^u) , 


and 

(ii)  Q(x,  y,  u,  v)  =  y^C'^v  +  (x  -  BC_‘*y)^(A-BC_‘*'Br^)~“'"(u  -  BC-^  v) . 
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PROOF:  The  positive  definiteness  of  A  and  C  is  obvious. 
Writing 


'A  B‘ 

'a  o' 

'I  A-1B 

T 

T 

T  -1 

Lb  cj 

Lb  ij 

Lo  C-B  A  Bj 

T  —X 

and  taking  determinants  of  both  sides,  we  see  that  |C-B  A  I  >  0, 

T  -1  -1 

which  shows  that  (C-B  A  B)  exists.  Inverting  both  sides  of  the 
above  decomposition  then  gives 


'A  B' 

-l 

_1  T  1  — 1  . 

I  -A  B(C-B  A  B) 

'  A-1 

0  ' 

T 

Lb  cj 

.  o  (c-bta-1b)-1. 

.-bta-1 

I . 

If  Q(x,  y,  u,  v,)  is  computed  using  the  representation  of  M-'*' 

given  by  the  above  equation,  (i)  is  obtained.  Setting  x  =  u  =  0, 

T  -1 

y  =  v,  the  positive  definiteness  of  C-B  A  B  is  apparent.  The 
proof  of  the  remaining  statements  is  similar. 

A/ 

The  maximum  likelihood  estimate  of  a  given  Z^is  simply  the 
value  of  a  which  minimizes  the  quadratic  form 


D*.§] 

>2  ' 

1 — 1 

1 

'a  ' 

“n  - 
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which  occurs  in  the  exponential  factor  of  the  joint  Gaussian  density 


of  a  and  Z^.  Using  expansion  (ii)  of  the  above  lemma,  this 
quadratic  form  may  be  written  as 


(3-(Eg^TUN~%)2 

(  <r2-(Ea^)TuN-1(E^^)) 


(V.7) 


From  this  last  expression,  we  see  that  the  maximum  likelihood  estimate 
of  a  given  2^j  under  the  Gaussian  assumptions  is  just  the  same  as 
the  general  linear,  minimum  variance  estimate  (V.5),  and  the  estimation 
error  also  has  the  same  variance  (V.6). 

Returning  now  to  the  wide-sense,  minimum- variance  point  of  view, 
we  note  that  not  only  for  9f,  but  in  fact  for  any  element  K  ,  the 
projection 


-  o*V)  y% 


(V.8) 


is  the  optimum  linear  estimate  of  K  given  ^  in  the  sense  of  the 
norm  of  W  ,  i.e.,  in  the  minimum-variance  sense.  If  y  is  a  random 
vector  such  that  the  components  of  the  associated  vector  y  are 
elements  of  X  ,  we  will  denote  by  that  vector  whose  components 
are  the  projections  on  (N)  of  the  components  of  y.  Hence,  we 
can  write 


y  =  o^T>  y%- 


(V.9) 
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The  covariance  matrix  of  the  estimation  error  will  then  be 
given  by 

[tf1  -  It/,  i  -  it/)] 

“  cov  (2)  -  (EJ^T)  Ujj-1  (e2n2T),  (V.10) 

where  and  $  denote  the  i—  and  j—  components  of  the  vector  y. 

The  minimum^ variance,  linear  estimate  of  the  original  vector  y  is 
then  y  +  K^y,  and  the  covariance  matrix  of  the  error  involved  in 
this  estimate  is  also  given  by  (V.10). 

We  will  now  re-write  the  system  (V.la)  in  the  augmented 

form 


y(k+l)  =  <5  (k+1,  k)  y(k)  +  u(k)  +  w(k) , 

J 


(V.ll) 


where 


Z(k)  = 


and 


$y(kfl,k) 


$(lcfl,k) 

0 


l(k) 

1 
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If  M(k)  is  defined  as  the  m  x  (n+l)  matrix  [H(k),  0  ]  ,  we  then 
have 


z(k)  =  M(k)  y(k)  +  v(k) 

(V.12) 

instead  of  equation  (V.2).  Together  with 

have  the  associated  equations 

(V.ll)  and  (v.12),  we 

y(k-KL)  =  (|)  (k+l,k)y(k)  +  u(k), 

V 

(V.lla) 

yCk+l)  =  (|)  (k+l,k)y(k)  +  w(k) 

«✓ 

(V.llb) 

z(k)  =  M(k)y(k),  and 

(V.l2a) 

z(k)  =  M(k)y(k)  +  v(k). 

(V.l2b) 

The  optimal  estimate  of  y(k)  given  JL  =  will  be  denoted 

by  y(k/N);  hence, 


y(k/N)  =  y(k)  +  K^k), 


(V.13) 


since  the  components  of  y(k)  are  obviously  elements  of  % .  The 
estimation  error  associated  with  this  estimate  will  be  designated 
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as  2(k/N)  ,  so  that 


£(k/N)  =  j;(k)  -  $(k/N) 

=  y(k)  -  K^(k).  (V.14) 

setting  P  (k/N)  =  cov  (^(k/N)),  it  follows  from  equation  (V .10) 

77 

that 


Pjjfk/h)  =  cov  (g(k))-  (is;(k)4T)  O^1  (E^2T(k)>.  (V.15) 

Since  w(k)  and  Z,  are  orthogonal  for  each  k,  we  compute 


l(k/k-l) 


$y(k,k-l)£(k-l)  +  u(k-l)  +  K^_^(  (£y(k,k-l)y(k-l)  +  w(k-l)) 

$  (k,k-l)£(k-lA-D  +  u(k-l).  (V.16) 

V 


Subtracting  (V.16)  from  (V.ll)  (after  replacing  k  by  k-1  in  (V.ll))> 
and  evaluating  the  covariance  matrix  of  the  resulting  expression, 
produces 


P  (k/k-l)  =  <J>  (k,k-l)P  (k-l/k-1) 
yy  *y  ’  y y 


$yT(k,k-l)  + 


Q(k-l) 


0 


(V.17) 

0 

0 
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Applying  expansion  (i)  of  the  above  lemma,  we  now  compute 


K^(k)  =  (E2(k)^T) 


=  D#(k)4:r 


EyCk)!1^)] 

'Uk-1 

E4_if(k) ' 

-1 

z,  . 

„  EfCk)^^ 

E|(k)|T(k)_ 

J(k)_ 

=  Kk_xy(k)  +  By(k)  (z(k)  -  ^Kk)), 


(V.18) 


where 


By(k)  -[E2(k)2T(k)  -  (E2(k)^L)  U-^CE^/Ck))] 

*[E^(k)zT(k)  -  (E^k)^)  U"^  (E^z^k))]-1. 


Since  v(k)  is  orthogonal  to  both  y(k)  and  Zk  we  can  substitute 
the  right  side  of  (V.12b)  for  2f(k)  in  the  latter  expression,  and 
obtain 


Bv(k)  -  P  (k/k-l)MT(k)  [M(k)P  (k/k-l)MT(k)  +  R(k)  ]  1 


yy 


yy 


(V  .19) 


The  orthogonality  of  v(k)  and  also  implies  that  Ft  ,z(k) 

MCk)!^  ny(k);  hence,  the  combination  of  (V.13)  and  (V.18) 
produces 

y(k/k)  =  y(k)  +  K^J^k)  +  By(k)  [  z(k)  -  M(k)[y(k)  +  Kk_1^( k ) J ] 

=  y(k/k-l)  +  By(k)  [  z(k)  -  M(k)y(k/k-l)  ]  .  (V.20) 
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From  (V.Hb),  (V.l2b),  (V.14)  and  (V.18),  we  find 


J(kA)  ■  z(k)  -  ^^(k)  -  By.(k)M(k)  [y(k)  -  K^^k)]  -  By(k)v(k) 

=  (I-B  (k)M(k));jr(k/k-l)  -  B  (k)v(k).  (V.21) 

J  J 

Computing  the  covariance  of  the  latter  expression  gives 

P  ( k/k )  =  ( I-B  ( k )M( k ) )  P  ( k/k-1 ) .  ( V . 22) 

V  J  J  V  J 

Equations  (V.16),  (V.17),  (V.19),  (V.20),  (V.21)  and  (V.22)  are 
equivalent  to  those  originally  derived  by  Kalman^  for  the  sequential 
estimation  of  the  state  vector  of  the  system  (V.ll)  based  on 
observations  of  the  form  (V.12)  (c.f.  equations  (3-5),  (3*6),  (3.14), 
(3.15),  (3*16)  and  (3.17)  of  reference  ll).  In  order  to  start  the 
computation,  it  is  clear  that  we  should  set 


p  (0/0)  = 
yy 

'po  0 

,  and  y(0/0)  = 

'x(0  )* 

1 - 

O 

t _ 

.  «  . 

(V.23) 


As  the  computation  proceeds,  the  optimal  estimate  of  the  parameter 
a  based  on  N  observation  points,  which  is  given  by  equation  (  .5), 
is  obviously  the  same  as  the  last  component  of  the  vector  y(N/N). 
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APPENDIX  VI 


THE  PSEUDO- INVERSE  AND  DATA  EDITING 


A.  INTRODUCTION 

The  optimal  noise-free  solution  for  the  Kalman  filter  derives  from 
finding  the  gain  B  (n)  that  satisfies  the  equation  (4*39): 

y 

B  [h  P  (n/n-l)  H  T 1  -  P  (n/n-1)  H  1  =  0,  (VI.l) 

y  L  y  yy  y  J  yy  y 


where  the  argument  (n)  is  implicit  in  the  B  and  H  terms.  In  order 
for  the  processing  of  perfect  observations  to  be  meaningful,  we  can  use 
a  maximum  of  only  three  position  measurements  and  three  velocity  measure¬ 
ments  at  any  time  t  ,  for  otherwise  we  would  have  some  wholly  redundant 


•  •  • 


equations  in  the  unknowns  x,  y,  z,  x,  y,  z  without  enhancing  our  knowledge 


of  the  mass  parameters  u^,  ur  ,  u^.  As  long  as  we  are  dealing  with  a  single 

sensor  at  time  t  ,  as  long  as  t  is  not  identical  to  t  ,  for  all  n  (no 
nJ  n  n-1 

matter  how  close  they  may  get),  and  as  long  as  the  trajectory  is  randomly 


perturbed  (with  geopotential  and  drag  uncertainties),  then  H  P  ( n/n-1 )H 

is  theoretically  nonsingular  and  the  solution  (4.40) 


B  (n)  =  P  ( n/n-1  )H  T[hP  ( n/n-1  )H  T1  ~1  (VI. 2) 

y  yy  y  l  y  yy  '  y  J 

is  theoretically  possible. 

In  a  practical  sense,  however,  the  inversion  called  out  in  (VI. 2) 

often  meets  with  severe  numerical  difficulties.  The  first  of  these 

arises  because  we  may  not  be  able  to  keep  t  distinct  from  t  .  for 

y  n  n-1 

all  n. 
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1.  Multiple-Sensors,  Simultaneous  Observations 


Two  or  more  sensors,  each  perfect,  may  be  able  to  see  the  satellite 


at  the  same  instant  of  time  t  .  For  any  of  the  sensors  of  concern  to  us 

n 

T 

in  the  present  work,  H  P  (n/n-l)H  would  then  be  singular. 

y  77  7 

EXAMPLE  1. 


Suppose  two  Baker-Nunn  cameras  sight  a  satellite  at  the  same  time 

t  .  Each  of  the  cameras  has  an  observation  matrix 
n 


da/d£ 

hll  h12  °  OOOOOOO' 

d£/dy 

,h2l  h22h23°00000° 

where  a  = 

s  = 


hll  = 


hl2 


h2l 


h0„  = 

22 


h23 


where  (x  , 
s’ 


right  ascension, 
declination, 

yy,  _ 

(sec2a)(x-x  )2 
s 

1 _ 

(sec2a)(x-x  )2 
s 

-(z-zs)(x-xs) 

RJ  cos  £ 

-(z-z g)(y-ys) 

O  * 

cos  £ 

R2  -  (z-z  )2 

:  _  3 

R^  cos  & 

y  ,  zg)  are  the  coordinates  of  the  camera  in  question, 


R  is 
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the  range  from  the  camera  to  the  satellite,  and  the  h.  ^  are  evaluated 

around  nominal  satellite  values  at  t  . 

n 

Vie  can  let  t  be  itself  for  camera  1,  and  let  it  be  t  for  camera 
n  ’  n+1 

To  make  the  computations  easier,  we  will  take  both  cameras  at  the  same 

latitude  such  that  z  =  z  ,  and  normalize  where  we  can  to  obtain 

s  ’ 


2. 


Camera  1: 


Hyi(n) 


110000000 

001000000 


Camera  2:  H  (n+1)  =  H  (n) 


abOOOOOOO 

OOcOOOOOO 


We  will  assume,  without  any  real  restriction,  that 


P  (n/n-l)  =  O'2!. 

yy 

Then  we  can  show  that 


(VI. 3) 


P  (n/n)  =  (<r2/2) 

*/  V 


1 

-1 


-  1  ! 


H 


10 - - 

I  2 
I  2 
I  2 


(VT.4) 


(This  is  also  P  (n+l/n)  since  the  data  from  camera  2  occurs  at  t  ,,  =  t  .) 

yy  n+1  n 

Mote  that  the  two  angles  from  camera  1  removed  two  degrees  of  freedom 
(reduced  the  rank  by  2)  between  (VI.3)  and  (VI. 4).  This  occurred  where 
it  should,  in  the  3x3  submatrix  in  the  upper  left-hand  corner,  which 
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corresponds  to  position  uncertainties.  To  remove  the  final  uncertainty 
in  that  3x3  submatrix,  we  should  use  only  one  more  angle  measured  at  that 
time:  the  right  ascension  reading  from  camera  2.  Without  a  special 
procedure,  however,  we  try  to  use  both  new  angles: 


Hy2(n+l)Pyy(n+l/n)Hy2T(n+l)  =  H^n+l^n/n^Cn+l) 


=  (<£/2) 


(a-b)' 


(VI. 5) 


and  we  obtain  an  expression  which  cannot  be  inverted  fcr  (VI. 2). 

2.  Round-Off 

Even  when  the  matrix  is  theoretically  well-behaved,  the  finite 
precision  of  the  computation  equipment  may  present  us  with  severe  numerical 
problems.  In  the  case  of  multiple  sensors,  two  observations  may  not  be 
truly  simultaneous ,  but  they  may  be  so  close  in  time  that  the  trajectory 
perturbations  have  had  essentially  no  effect  on  the  orbit:  i.e.,  the  obser¬ 
vations  are  treated  numerically  in  much  the  same  way  as  led  to  (VI.5).  Note 
that  the  perturbations  enter  the  covariance  computations  through  (^^(n-l) 
and  (n,  n-l),  both  of  which  are  integrals  of  finite  functions  over  the 

range  (t  ,  t  ,  ),*■  and  hence  vanish  as  t  1  ,. 

n  n-i  n  n-l 

Problems  can  still  arise  with  closely  spaced  observations  after  care 
has  been  taken  to  throw  away  redundant  readings,  such  as  the  declination 


*See  equations  (4.27)  and  (4-30). 
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measurement  in  the  above  example.  The  transition  equations  for 

the  stochastic  and  deterministic  estimates  of  drag  are  given  respectively 

as 

U-^(n/n-l)  =  u, (n-l/n-l)  e-^n~^n-l^"^^  , 

u,  (n/n-1)  =  u^n-l/n-l) .  (VI. 6) 

If  t  is  too  close  to  t  the  two  equations  are  numerically  difficult 
to  distinguish,  and  since  un  and  u  appear  in  the  equations  of  motion 
only  as  the  sum  (u^  +  u  ),  we  find  ourselves  faced  with  a  system  which  is 
essentially  Kalman-unobservable 1 
EXAMPLE  2. 

Consider  the  hypothetically  simplified  example 


where 


x(n+l)  =  j|  (n+1 ,n)  x(n), 
z.(n)  =  H(n)  x(n)  , 


$  (n+l,n)  =  |  = 


H(n)  =  H  = 


110  0 
0  111 
0  0  10 
0  0  0  1 

12  0  0 

3  4  0  0 


(VI. 7) 


Components  x^  and  x^,  here,  correspond  to  the  u^  and  u^  of  our  real  system 
under  the  condition  that  t^  and  t  are  so  close  together  that  the 
exponential  in  (VI. 6)  is  numerically  unity. 
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After  one  measurement,  say  at  tn, 


Cov 


*L 


=  0. 


We  will  assume  that  the  overall  covariance 

0 


P(n/n) 


\2j 


Then 


P(n+l/n)  =  $  P(n/n)  $  T  = 


and 


H  P(n+l/n)HT  =  U(^2+  \2) 


(T2+  or2 

5  4 

T2 


<r  2 

4 


3 

0 


<r2 

4 


(VI. 8) 


Again,  we  cannot  perform  the  inversion  required  for  (VI. 2).  What  we  should 
have  done  was  to  accept  only  one  of  the  two  observations  at  this  time 
and  estimate,  in  effect,  only  the  sum  of  and  x^,  rather  than  try  to 
estimate  them  independently.  Then,  when  we  obtain  some  better-spaced 
future  data  that  causes  their  behavior  to  separate,  estimate  them  individually 
B.  THE  PSEUDO- INVERSE 

Rather  than  attempt  to  throw  away  data  in  the  conscious  way  just  described 
we  can  achieve  the  same  effect  by  replacing  the  inverse  in  (VI. 2)  with  a 
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pseudo-inverse.  The  general  method  of  pseudo-inversion,  with  overtones 
for  least-squares  fitting,  is  described  in  Penrose, ^  and  Kalman"^ 
alludes  to  it  throughout  his  work  without,  however,  going  through  the 
mechanics  of  what  it  achieves. 

For  our  purposes,  since  we  require  it  only  for  the  observation- 

T 

covariance  matrix  I!  P  (n/n-l)H  ,  it  suffices  to  specialize  the  pseudo- 

y  yy  y 

inversion  to  a  symmetric  k  x  k  matrix,  say  A.  Denoting  the  eigenvalues 
of  A  as  k-^,...,X^,  we  can  always  transform  A  to  the  diagonal  form 


K 


D  = 


T 

=  5  AS, 


(VI. 9) 


by  taking  S  to  be  a  matrix  whose  k  columns  aro  the  k  distinct  eigen¬ 
vectors  of  A,  normalized  such  that 


T 

S  S 


T 


If  some  of  the  \fs  are  zero,  say  , , 

inverse  of  A  is  defined  as 


>  ,  then  the  pseudo- 


A*  -  (SDST)*  =  S 


(iAx) 


(tA2) 


(l/xp 


o 


o 
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That  is,  the  part  of  A  that  has  a  normal  inverse  is  inverted  normally; 
the  part  that  has  zero  eigenvalues  is  pseudo-inverted  to  have  zero 
eigenvalues . 

APPLICATION  TO  EXAMPLE  1. 

Consider  the  solution  with  the  declination  measurement  thrown  away 
by  conscious  editing.  We  use  (  )  to  denote  the  matrices  so  obtained 
for  camera  2: 

z '  =  a, 

y 

H-  -  [a  b  0  0  0  0  O  0  O] 


Then,  from  (VI. 2),  (VI.4),  and  (VI.5) 

B'(n+1)  =  P  (n/n)H'^(n+l)[H'  (n+l)P  (n/n)H^(n+l)]  ", 


-1 


=  (<r2/  2) 


a-b 

'  1  ' 

a+b 

a+b 

,  -1 

1 

0 

# 

• 

• 

l(<r2/2)(a-b)2]  = 

a-b 

0 

• 

# 

0 

• 

• 

0 

(vi.10) 


Now  we  apply  the  pseudo-inverse,  instead,  to  (VI. 5)  in  the  automatic 
way  described: 

Vn+1)  -  Vn/n)Vn+1)  [  Hy2(n+1) Vn/n)Hy2(n+1)]# 
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=  (<r2/ 2) 


a-b 

a+b 

0 

♦ 

♦ 

0 


((T2/  2)  (a-b) 2  0 


0 


1 

a+b 


a-b 

0 

# 

0 


(VI. 11) 

Relation  (VI.10)  directs  that  the  first  two  position  extrapolations  be 
corrected  by  multiplying  the  right-ascension  residuals  by  l/(a+b)  and 
l/(a-b),  respectively,  and  all  others  be  left  uncorrected;  (VI. 11 )  directs 
the  same  thing.  Hence  the  results  are  the  same. 

APPLICATION  TO  mHPLL  2. 

We  accept  only  the  z-,  measurement  in  the  simplified  example  specified 
by  (VI. 7).  Then 

I 

Z  =  Z-, 


n'  =  [1  200] 


and 


B  (n+1)  =  2 


0 

0 

9  9 

cC,  +<r  ^ 

CT2  +77  * 

3  4 

2  2  1  1 

3  4 

o;2 

(4^  4  4^  )  “ 

3  4 

r  2 

3 

<r  2 

<r.2 

4 

4 

(VI.12) 
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The  in-step  covariance  is 

P' (n+l/n+1 )  =  [l-B'(n+l)H'JP(n+l/n), 


=  _2 


<r32+^2 


0  0  0  0 

0  0  0  0 

0  0  1-1 

0  0-11 


(VI. 13) 

The  pseudo-inverse  approach  employs  the  original  unprimed  matrices 
from  Example  2.  First  we  must  find  the  eigenvalues  of  (VI. 8): 

X^  =  5,  X^  =  0.  We  then  determine  the  matrix  of  eigenvectors  to  be 


-2 


s  =  sT » <i//r  ) 


Equation  (VI. 8)  becomes 


H  P(n+l/n)n  =  20  *  +  <T^)  S 


so  that 


i 


0  0 


S, 


r  Tl* 

[H  P(n+l/n)H  J  = 


20(^2+<£2) 


100(0- 2+<r2) 
3  4 


0  0 
0  1 

1  2 

2  4 


(VI.14) 


Hence 


B(n+1)  = 


io(  tr32+cr42) 


0 

<tA<72 

J>  4 


0 

2(q-2+<T2) 
3  4 


2C^ 


K 


(vi.15) 
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Comparison  of  (VI. 15)  with  (VI. 12)  is  not  particularly  enlightening. 
Hence  we  compute  the  in-step  covariance  with  the  unprimed  values , 


P(n+l/n) 


0  0  0  0 

0  0  0  0 

0  0  1-1 

0  0-11 


and  find  that  it  agrees  with  (VI. 13),  thus  verifying  that  the  pseudo¬ 
inverse  automatically  performs  the  desired  editing* 


C.  NUMERICAL  IMPUEMENTAT ION 


Computer  round-off  usually  prevents  a  calculation  from  yielding  a  true 
zero  value  when  there  should  be  one  if  the  calculation  involves  more 
than  a  very  few  steps.  Hence  the  examination  of  D  for  computed  singularities 
will  not  always  be  sufficient. 

There  are  two  cases  to  consider.  The  first  and  most  common  arises 

2 

when  some  diagonal  entry  <T  in  should  be  zero,  but  is  not.  The  entries 

2 

in  the  row  and  column  that  contain  cr  should  also  be  zero  if  P  is  to 

i  yy 

have  the  non-negative  definite  property  that  a  covariance  matrix  must  have. 

If  these  row  and  column  entries  are  wholly  due  to  round-off,  it  is  very 
unlikely  that  P  will  be  non-negative  definite  (have  only  non-negative 

yy 

eigenvalues).  Hence,  one  approach  is  to  positive-semi-definitize  the  P 

yy 

matrices . 


TEST  1. 


Check  the  diagonal  entries  in  P 

yy 


If  there  are  any  negative  entries, 
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say  »  set  them  and  the  ith  column  and  row  to  zero.  Set  also  the 
appropriate  rows  in  Py,j  and  P^  .  to  zero.  Perform  the  check  for  both 
the  arguments  (n/n)  and  (n/n-l). 

TEST  2. 

Compute  the  diagonal  matrix  of  eigenvalues 
D 

for  arguments  (n/n)  and  (n/n-l).  If  any  eigenvalue  is  negative,  set  it 


T 

S  P  S 
77 


to  zero  and  reconstruct 

P  '  =  SD'ST, 

yy 

where  the  primes  denote  the  matrices  after  the  negative  eigenvalues  are 
made  zero.  It  is  assumed  that  any  negative  eigenvalue  is  small  enough 
in  magnitude  that  setting  it  to  zero  has  negligible  effect  on  S. 

TEST  3. 

For  some  input  parameter  C,  set  the  ith  row  and  column  in  P  (n/n) 

77 

to  zero  if  tT^2  C. 

TEST  4. 

Having  obtained  a  data  value  z. (n),  and  therefrom  a  residual 


ei  =  zi(n)  -  hi[x(n/n-l)j  , 


set  the  ith  element  <f  .  in  H  (n)P  (n/n-l)H  (n)  to  zero  if 
—  ei  7  77  7 

e2  >  K,  <T  2 

1  1  ei 

for  some  input  parameter  K^.  However,  set  the  residual  e^  to  zero  if 

2  — 2  -K  5T2. 


\  <T  r  >  e  r  > 


ei  -  i 


2  ei 
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for  some  input  parameter  K,  (<  K^)  and  proceed  with  the  covariance 
computations  as  if  the  data  was  used. 

The  second  case  arises  when  the  P  matrices  are  numerically  positive 

yy 

semi-definite ,  but  some  small  should  be  zero.  There  are  two  possible 

approaches.  One  is  to  diagonalize  to  the  D  matrix  and  set  any  X^  to 

zero  that  is  Mvery  much"  smaller  than  any  other  There  is  a  difficulty 

in  assigning  a  number  to  the  "very  much",  especially  since  some  of  the 

variances  have  the  units  of  distance,  some  have  the  units  of  velocity,  and 
some  have  units  involving  area-to-mass  ratio.  In  general,  ha*/ ever,  it 
is  better  to  discard  good  data  (by  setting  a  X^  to  zero  that  should  not 
be)  rather  than  to  include  bad,  or  meaningless  data  in  the  smoothing. 

The  computed  covariances  will  simply  be  a  little  larger  than  they  should 
be . 

A  second  approach,  which  also  covers  some  of  the  problems  encountered 

earlier,  is  to  include  a  round-off  "noise"  matrix  R  in  the  residual 
’  rr 

covariance 

tHy(n)Pyy(n/n-l)Hy(n)  + 

to  reflect  the  fact  that  we  do  not  have  an  infinitely  precise  processing 
chain,  even  when  the  first  unit  in  the  chain  (the  sensor)  is  perfect. 


-  166  - 


UNCLASSIFIED _ 

_ Security  Classification 


DOCUMENT  CONTROL  DATA  -  R  &  D 


(Security  classification  of  tltla,  body  of  abstract  and  indexing  annotation  must  be  entered  when  tha  ovaratl  report  te  classified) 


1 .  ORIGIN*  TING  AC  Tt  VI  TY  (Corporate  author) 

Westing  house  Defense  and  Space  Center 

Surface  Division 

Baltimore,  Maryland 

2a.  REPORT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 

2b.  GROUP 

N/A 

3.  REPORT  TITLE 

MODELS  FOR  ANALYSIS  OF  THE  CAPABILITIES  OF  G  ROUND  BASED  SENSORS 

IN  DETERMINING  THE  MASS  OF  ORBITING  BODIES 

4  DESCRIPTIVE  NOTES  (Type  of  report  and  Indus iva  datas) 

None 

5 ■  AUTHOR (5)  (First  name,  mtddia  initial ,  last  name) 

Walter  J.  Culver,  et  al 

8  REPORT  DATE 

23  June  1967 

7a.  TOTAL  NO.  OF  PAGES  7b.  NO.  OF  REFS 

166  50 

8a.  CONTRACT  OR  GRANT  NO 

F19628-67-C0041 

b.  P  ROJ  EC  T  NO. 

C. 

d. 

9a.  ORIGINATOR'S  REPORT  NUMBER(S) 

ESD-TR-68-157,  Vol.  1 

9b.  OTHER  REPORT  NOtsI  (Any  other  numbers  that  may  be  assigned 
this  report) 

None 

10.  DISTRIBUTION  STATEMENT 

This  document  has  been  approved  for  public  release  and  sale;  its  distribution  is  unlimited. 

11  SUPPLEMENTARY  NOTES 

12.  SPONSORING  MILITARY  ACTIVITY 

Space  Defense  Systems  Program  Office,  Electronic 
Systems  Division,  Air  Force  Systems  Command, 
USAF,  LG  Hanscom  Fid,  Bedford,  Mass.  01730 

1  3  ABSTRAC T 


This  report  contains  a  description  of  the  models  to  be  used  in  analyzing  the  capabilities  of 
ground-based  sensors  h  determining  the  mass  of  orbiting  bodies,  model  coefficients,  and  the 
justification  for  their  selection.  Relations  are  derived  for  computing  sensitivity  coefficients 
and  their  coupling  to  mass  variance. 


DD  ,Fr,,1473  UNCLASSIFIED 


Security  Classification 


UNCLASSIFIED 


Security  Classification 


1  4. 

LINK  A 

LINK  B 

LINK  C 

KEY  WO  RDS 

ROLE 

W  T 

ROLE 

W  T 

ROLE 

W  T 

RADAR 

MASS 

ERRORS 

ANALYSIS 

SENSOR 

ORBIT 

SATELLITE 

UNCLASSIFIED 


Security  Classification 


